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A  digital-computer  program  has  been  written  for  calculating  the  viscous 
drag  of  streamlined  bodies  of  revolution  In  constant-density  axial  flow.  The 
Integral  approach  adopted  Incorporates  recently  Improved  methods  for  pre¬ 
dicting  the  transition  point  and  for  calculating  the  turbulent  boundary 
layer.  The  Inputs  to  the  computer  program  are  the  body  geometry,  the 
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ABSTRACT 


A  digital-computer  program  has  been  written  for  calcu¬ 
lating  the  viscous  drag  of  streamlined  bodies  of  revolution 
in  constant-density  axial  flow.  The  Integral  approach  adopted 
incorporates  recently  Improved  methods  for  predicting  the 
transition  point  and  for  calculating  the  turbulent  boundary 
layer.  The  inputs  to  the  computer  program  are  the  body 
geometry,  the  associated  invlscid  pressure  distribution,  and 
the  body-length  Reynolds  number.  Agreement  of  calculated  and 
measured  drag  coefficients  is  good,  particularly  in  cases 
where  the  transition  point  is  predicted  accurately. 

ADMINISTRATIVE  INFORMATION 

Tlie  work  reported  here  was  supported  by  the  in-house  Independent 
research  program  of  the  David  W.  Taylor  Naval  Ship  Research  and  Develop¬ 
ment  Center  (DTNSRDC)  and  funded  under  Task  Area  ZR-023-0101,  Work  Unit 
1-1541-002. 


INTRODUCTION 

The  accurate  prediction  of  the  viscous  drag  of  streamlined  bodies 
of  revolution  in  steady  axial  motion  in  a  constant-density  fluid  is  a 
basic  hydrodynamic  problem  of  great  importance  to  designers.  A  stream¬ 
lined  body  is  defined  here  as  one  on  which  there  is  no  significant  flow 
separation;  thus  such  configurations  as  those  with  very  blunt  sterns  or 
noses  are  excluded.  In  1953  Granville^  reported  a  method  of  calculating 
viscous  drag  based  on  momentum- integral  formulations  of  the  laminar-  and 
turbulent-boundary-layer  equations.  The  method  of  predicting  the  transi¬ 
tion  point  was  based  on  a  correlation  of  two-dimensional  airfoil  data  in 
terms  of  the  difference  between  the  momentum-thickness  Reynolds  numbers 
at  the  transition  and  neutral-stability  points  and  the  averaged 


Granville,  P.  S.,  "The  Calculation  of  the  Viscous  Drag  of  Bodies  of 
Revolution,"  David  Taylor  Model  Basin  Report  849  (Jul  1953).  A  complete 
listing  of  references  is  given  on  page  141. 
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pressure-gradient  parameters  between  these  two  points.  The  turbulent- 
boundary-layer  calculations  were  based  on  power-law  relations  which  result 
In  simple  quadratures.  The  calculations  were  separated  Into  two  parts, 
according  to  whether  the  boundary  layer  was  thick  or  thin  relative  to  the 
local  body  radius. 

The  updated  method  of  drag  prediction  described  herein  follows  that 

of  Granville^  In  outline.  The  laminar-boundary-layer  calculations  are 

unchanged.  The  transition  point  Is  found  from  a  new  correlation  of 

transition  data  for  bodies  of  revolution  in  terms  of  the  difference 

between  the  momentum-thickness  Reynolds  numbers  at  the  transition  and 

neutral-stability  points  and  the  rate  of  change  of  body  shape.  This  new 

2 

correlation  was  obtained  by  Granville  and  applies  to  curved  bodies;  as 

an  alternative,  the  method  of  transition  prediction  obtained  earlier  by 
1  3 

Granville  and  applied  by  Smith  can  be  used.  The  turbulent-boundary- 

layer  calculations  are  made  from  two  "integral"  differential  equations, 

a  momentum  equation  and  a  shape-parameter  equation  based  on  entrainment. 

A 

The  new  method  is  developed  by  Granville  and  incorporates  velocity- 
similarity  laws.  It  applies  to  either  thick  or  thin  boundary  layers  and 
so  it  is  not  necessary  to  give  separate  consideration  to  the  region  near 
the  tail  where  the  boundary  layer  is  thick. 

A  digital-computer  program  has  been  written  to  execute  the  calcula¬ 
tions  and  compute  the  drag.  A  preliminary  program  is  used  to  compute  the 
body  shape  from  a  specified  polynomial  description;  this  is  used  as 
input  to  the  Douglas-Neumann  potential-flow  program  which  calculates  the 
inviscid  pressure  distribution  for  use  in  the  boundary-layer  calculations. 


2 

Granville,  P.  S. ,  "The  Prediction  of  Transition  from  Laminar  to 
Turbulent  Flow  in  Boundary  Layers  on  Bodies  of  Revolution,"  NSRDC  Report 
3900  (Sep  1974);  also  presented  at  the  Tenth  Office  of  Naval  Research 
Symposium  on  Naval  Hydrodynamics,  Massachusetts  Institute  of  Technology, 
Cambridge,  Mass.  (Jun  1974)  and  will  appear  ia  the  proceedings  of  this 
symposium. 

3  I 
Smith,  A.  M.  0.,  Discussion  of  Granville  symposium  paper  ;  will 

appear  in  the  symposium  proceedings. 

4 

Granville,  P.  S.,  "Similarity-Law  Method  for  Thick  Axisymmetric 
Turbulent  Boundary  Layers  in  Pressure  Gradients,"  DTNSRDC  Report  4525 
in  preparation) . 
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Comparison  of  the  computed  results  with  towing-tank  measurements  shows 
(1)  that  the  accuracy  of  the  two  Granville  methods  for  predicting  transi¬ 
tion  is  roughly  equal  if  transition  occurs  on  the  forebody  and  (2)  that 
the  drag  is  predicted  accurately  if  the  transition  prediction  is  accurate. 
The  turbulent-boundary-layer  theory  on  which  the  calculations  are  based 
includes  the  Schoenherr  frictional  line;  the  program  can  be  forced  to 
reproduce  the  Schoenherr  line  as  its  predicted  drag  coefficient  by 
setting  the  pressure  gradient  equal  to  zero,  setting  the  body  radius 
equal  to  a  constant  sufficiently  large  value,  and  forcing  transition  at 
the  nose.  Granville  uses  the  Schoenherr  line  as  a  baseline  for  his 
method  because  of  its  classical  and  scientific  importance  in  turbulent¬ 
boundary-layer  theory.  In  the  figures  in  this  report  which  present 
calculations  done  by  this  theory,  the  Schoenherr  line  has  been  drawn  in 
for  comparison;  it  is  an  easy  matter  to  draw  in  other  lines  which  are  in 
widespread  use,  such  as  the  1957  International  Towing  Tank  Conference 
correlation  line  for  ship  models. 

A  method  for  predicting  the  viscous  drag  of  a  body  of  revolution 
has  also  been  reported  by  Nakayama  and  Patel. ^  It  is  similar  to  the 
method  of  Granville  reported  here  in  that  an  entrainment  equation  is 
used  (but  with  a  more  restricted  one-parameter  system)  and  in  that 
careful  consideration  is  given  to  the  region  near  the  tail  where  the 
boundary  layer  is  thick.  Four  alternative  methods  of  predicting  the 
transition  point  are  available,  one  of  which  is  that  of  Granville.^ 

Good  agreement  with  measured  results  is  reported. 

A  second  method,  reported  by  Cebeci,  Mosinskis,  and  Smith, ^  uses 
the  more  time-consuming  differential  formulation  of  the  boundary-layer 
equations,  with  an  eddy-viscosity  profile.  These  authors  provide  two 
alternative  methods  for  predicting  transition,  one  of  which  is  that  of 
Granville.^  Again,  good  agreement  with  measured  results  is  reported. 


Nakayama,  A.  and  V.  C.  Patel,  "Calculation  on  the  Viscous  Resistance 
of  a  Body  of  Revolution,"  Journal  of  Hydronautics ,  Vol.  8,  No.  A, 
pp.  154-162  (Oct  1974). 

^Cebeci,  T.  et  al.,  "Calculation  of  Viscous  Drag  in  Incompressible 
Flows,"  Journal  of  Aircraft,  Vol.  9,  No.  10,  pp.  691-692  (Oct  1972). 


3 


Parsons,  Goodson,  and  Goldschmled ^  used  the  method  of  Cebecl  et  al.^ 
together  with  an  optimum-search  strategy  to  find  the  body  of  revolution 
which  has  minimum  drag  for  specified  speed  and  enclosed  volume.  Their 
approach  is  to  find  a  shape  which  has  transition  as  far  downstream  as 
possible  at  the  specified  speed.  Good  agreement  with  measured  results 
is  reported. 

g 

Following  recent  German  methods  described  by  Lugt  and  Oh,  Oh  and 
Reingruber*  have  used  an  Integral  formulation  of  the  boundary-layer 
equations  and  a  transition  criterion  different  from  those  reported  here 
to  compute  the  viscous  drag  on  a  number  of  bodies  of  revolution.  The 
only  comparison  with  measured  resucs  provided  is  for  the  transition 
point  on  one  body,  and  fairly  good  agreement  was  obtained  in  that 
instance. 

DESCRIPTION  OF  DIGITAL-COMPUTER  PROGRAM 

The  digital-computer  program  which  is  the  subject  of  this  section 

of  the  report  predicts  the  viscous  drag  on  a  streamlined  body  of  revolu- 

1  4 

tion  according  to  the  boundary-layer  theory  developed  by  Granville.  ’ 

The  program  actually  consists  of  three  distinct  programs;  if  desired, 
the  execution  of  one  of  these  can  be  partially  repeated  and  the  execution 
of  a  second  can  be  completely  repeated  for  each  body  and  length  Reynolds 
number  R^^.  These  three  programs  are  as  follows: 

DPIN  This  program  computes  the  shape  (offsets,  slopes,  etc.) 

of  the  body  of  revolution  on  the  basis  of  geometric 
parameters  which  are  input  in  accordance  with  polynomials 
developed  by  Granville.  At  present  three  alternative 
DPIN's  are  available  (DPINl,  UPIN2,  and  DPIN3)  for  three 

^Parsons,  J.  S.  et  al.,  "Shaping  of  Axisymmetric  Bodies  of  Minimum 
Drag  in  Incompressible  Flows,"  Journal  of  Hydronautics,  Vol.  8,  No.  3, 
pp.  100-107  (Jul  1974). 

g 

Lugt,  H.  J.  and  S.  K.  Oh,  "Boundary-Layer  Suction  with  Slots  on  Axi¬ 
symmetric  Bodies,"  NSRDC  Report  4038  (Nov  1972). 

*As  reported  Informally  in  NSRDC  Technical  Note  CMD-30-73  dated 
September  1975. 


DA  50 


DPOUT 


d  i  t  f ortMit  i  larses  of  bodies;  the  user  may  create  addi¬ 
tional  alternatives  for  other  class  of  bodies.  In  this 
report,  when  DPIN  is  mc'ntloned,  It  Is  to  be  understood 
that  the  particular  DPIN  is  meant  which  describes  the 
body  being  treated.  Appendixes  A-C  provide  descriptions 
of  the  respective  DI’lN’s. 

This  is  the  Doug las-Neumann  program  developed  by  Smith, 
Hess,  and  their  associates  tor  computing  inviscid  sub¬ 
sonic  flow  about  a  body  of  revolution.  It  represents  the 
body  by  a  series  of  frustums  of  cones,  with  axes  along 
the  axis  of  symmetry,  and  it  assumes  a  constant  (unknown) 
hydrodynamic  source  strength  on  each  frustum.  The 
boundary-value  problem  for  the  Laplace  equation  with  the 
boundary  condition  of  zero  normal  velocity  on  the  body 
leads  to  a  Fredholm  integral  equation  of  the  second  kind. 
This  is  solved  as  a  set  of  linear  algebraic  equations  nv 
using  (in  the  option  selected  here)  Seidel  iteration. 

The  velocity  potential,  velocity  components,  and  pressure 
anywhere  in  the  flow  field  can  be  found  relatively  easily 
after  the  distribution  of  source  strength  on  the  bodv  is 
determined. 

This  program  computes  the  boundary  layer  on  the  body  and 
the  viscous  drag,  using  the  shape  computed  by  DPIN  and 
the  inviscid  velocity  distribution  computed  by  DA50.  It 
also  computes  the  boundary- layer  displacement  thickness 
and  relates  an  increment  in  hydrodynamic  source  strength 
to  the  (inviscid)  velocity  and  the  displacement  thickness 
Then,  if  desired,  the  program  can  go  back  to  the  end  of 
DA50,  add  in  the  increment  in  source  strength,  recompute 
the  velocity  distribution,  return  to  DPOUT,  and  recompute 
the  boundary  layer  to  give  a  more  accurate  prediction. 
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In  principle  this  iteration  cruld  be  repeated  until  some 
convergence  criterion  is  satlsiled;  however,  in  Its 
present  form,  the  program  recomputes  the  boundary  layer 
only  once.  The  advantage  of  this  method  for  incorpor¬ 
ating  the  effect  of  displacement  thickness  is  that  the 
return  to  DA50  is  at  a  point  following  the  time-consuming 
Seidel  iteration.  Thus  this  iteration  need  be  executed 
only  once  for  a  given  body  of  revolution  even  though  the 
boundary  layer  is  computed  at  several  each  of  which 

would  correspond  to  a  different  displacement  thickness. 
Appendix  D  provides  descriptions  of  DA50  and  DPOUT. 

Figure  1  illustrates  a  representative  body  of  revolution  with  its 
boundary  layer  in  the  coordinate  system  used  in  the  calculations. 

CALCULATION  OF  LAMINAR  BOUNDARY  LAYEF 

The  method  by  which  the  laminar  boundary  layer  is  calculated  has 

1 

been  described  by  Granville.  For  the  case  of  zero  pressure  gradient, 

the  integral  method  used  reduces  to  the  Blasius  flat-plate  relation. 

The  possibility  of  laminar  separation  is  considered  according  to  the 

9 

criterion  of  Thwaites.  The  increment  in  source  strength  used  to  repre¬ 
sent  displacement  thickness  is  calculated  according  to  Lighthill.^^ 

The  method  begins  by  calculating  the  pressure-gradient  parameter 

2 

6**  dU  dU 

- ;  where  9,  -r-,  and  v  are  the  dimensional  momentum  thickness,  velo- 

V  ds’  ’  ds’ 

city  gradient,  and  kinematic  viscosity  at  points  along  the  laminar 

2 

0^  dU 

boundary  layer.  The  FORTRAN  name  for  ”  ^  DPOUT  is  PGF.  Equa¬ 
tion  (52)  in  Granville^  gives 


9 

Thwaites,  B.,  "Approximate  Calculation  of  the  Laminar  Boundary 
Layer,"  Aeronautical  Quarterly,  Vol.  I,  Part  III,  pp.  245-280  (Nov  1949). 

^^Lighthill,  M.  J.,  "On  Displacement  Thickness,"  Journal  of  Fluid 
Mechanics,  Vol.  4,  Part  4,  pp.  383-392  (Aug  1958). 


6 


V 

k 


Fiaure  1  -  Body  of  Revolution  with  Regions  of  Potential,  Boundary-Layer. 

and  Wake  Flows  Indicated 

(Adapted  from  Figure  I  of  Granville’) 
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where  Rg  Is  the  Reynolds  number  based  on  U  and  0  and  R0/Rj^  is  given  by 
Equation  (51)  in  Granville^  as 


R 

^  _  4 _ B 

h.  ^  (U/UJ^  (y/L)^ 


(2) 


with 


Rg/R^  and  3  have  the  FORTRAN  names  RTH2RL  and  B  in  DPOUT.  According  tc 


Thwaites,^  laminar  separation  takes  place  at  the  point  where  first 

becomes  less  than  or  equal  to  -0.09. 

The  Blasius  flat-plate  shear-stress  coefficient  (named  CF  in 
DPOUT)  is  given  by 


X 

Cj  =  I*  y(x')  C^(x')  dx’ 

^  *0 


(4) 


Here  A^*"^  is  total  area,  C  (x)  *  and  R  is  the  Reynolds  number 

s 

based  on  and  distance  s  along  the  body  surface.  (The  axial  component 

of  is  cos  a  and  the  element  of  integration  is  ds  or  sec  a  dx.) 
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The  increment  In  source  strength  needed  Co  represent  the  effect  of 
displacement  area  A*,  defined  as 


A*  = 


(5) 


(6  Is  boundary  layer  thickness  and  n  is  coordinate  normal  to  surface) 
is,  according  to  Lighthlll,^^ 

sec  a  d 

a  =  Tp- —  T-  (U6*)  (6) 

2  TT  y  dx 

Here  0  is  the  source  strength  per  unit  area  normalized  such  that  the 
velocity  potential  due  to  a  fundamental  source  E  at  a  point  x^,  y^,  z^ 

Is 


4'(x,y,z) 


(x  -  x^)^  +  (y  -  y^)^  +  (z  - 


(7) 


For  the  laminar  boundary  layer,  the  power-law  profile  assumed  in 
Granville^  gives 


A*  =  3 


(8) 


where  is  the  momentum  area  defined  by 

o 

(6  and  n  arc  defined  under  Equation  (5)).  Equation  (35)  in  Granville^ 
is 

-  y  0  (10) 


9 


F' 


80 


y  ^ 

u 

R.  U 
0  00 


or 


LR^U 


(11) 


Rg  is  found  as  the  square  root  of  R^  times  the  right-hand  side  of  Equa¬ 
tion  (2),  and  Q/L^,  denoted  by  0ML2  in  DPOUT,  is  given  by  Equation  (11). 

2 

Equation  (8)  gives  A*/L  (ELSL2  in  DPOUT)  and  Equation  (6)  gives 


£_  sec  g  L  r  A*  *^^^^^0°^  U  d(A*/L^)  ' 

U  "  2  TT  y  ,2  d(x/L)  d(x/L) 

00  L  L  ®  J 


d(U/Uj 

d(x/L) 


has  been  calculated  near  the  beginning  of  DPOUT; 


tiating  subroutine  t)GT3  is  used^^  to  calculate  ^ 


—  is  denoted  by  DELSIG  in  DPOUT.  Calculation  of  DELSIG 
^00 


the  differen- 
The  quantity 
is  omitted 


when  execution  of  DPOUT  is  carried  through  the  second  time  for  a  given 
configuration  and  (through  use  of  the  control  variable  ICONTROL) 
since  there  is  no  further  need  for  it.  More  precise  values  of  the 
source  strength  are  not  needed  since  no  further  executions  of  DPOUT 
are  performed. 


^^"360  Scientific  Subroutine  Package,  Version  III,"  IBM  Reference 
Manual,  Serial  H-20-0205-3,  p.  319  (1968). 
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t 


Actual  calculation  of  —  by  using  Equation  (12)  Is  not  performed 

'^OO 


In  the  laminar-boundary-layer  part  of  DPOUT  unless  no  turbulent  boundary 
layer  Is  reached  on  the  body  due  to  very  small  values  of  Rj^.  Values  of 

A*  o 

~  are  calculated  by  using  Equations  (11)  and  (8),  and  ^  Is  found  from 

Equation  (12)  at  the  end  of  the  turbulent-boundary-layer  calculations. 

12 

Equation  (2.1)  In  Hess  and  Smith  shows  that  the  fundamental 
source  strength  in  the  Douglas-Neumann  programs  Is  the  negative  of 
Equation  (7).  Therefore,  the  calculated  DELSIG  must  be  subtracted 
from  rather  than  added  to  the  calculated  source  strength  in  LINK  6  of 
DA50. 

In  order  to  find  the  neutral-stability  point,  the  curve  of  Rg 
2 

9  ^  d  U  X 

versus  —  ^  (Figure  3  in  Granville  )  is  approximated  by  a  cubic  in 
2 

9^  dU 

—  — .  This  figure  and  the  numerical  approximation  to  it  are  shown  as 
Figure  2.  The  neutral-stability  point  is  found  by  comparing  Rg  with 
2 

0 **  dU 

the  function  of - r-  described  by  the  line  in  Figure  2  or  by  the  cubic  in 

V  ds 

the  program,  starting  at  the  nose  of  the  body.  The  neutral-stability 
point  is  assumed  to  be  the  point  nearest  the  nose  at  which  Rg  equals  or 


2 

0  **  dU 

exceeds  the  function  of - 

V  ds 

After  the  neutral  stability  point  is  found,  the  laminar-boundary- 

2 

9  ^  dU  2  2 

layer  calculations  are  continued.  ~  >  A*/L  ,  and  C  are 

V  as  D  t 


found  at  each  point.  A  function  of  body  geometry,  “  d  (x/L)  ’ 


12 

Hess,  J.  L.  and  A.  M.  0.  Smith,  "Calculation  of  Potential  Flow 
about  Arbitrary  Bodies,"  in  "Progress  in  Aeronautical  Sciences," 
Pergamon  Press,  Oxford  and  New  York,  Vol.  8,  (1966),  pp.  1-138. 
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Figure  2  -  Curve  for  Calculating  Ra'"'  as  a  Function  of  —  — 

ds 

and  Numerical  Approximation  to  This  Curve 


12 


denoted  by  TP,  is  also  calculated  at  each  point  dovmstreani  of  the 

neutral-stability  point.  (D^  Is  maximum  body  diameter.)  TP  will  be 

needed  if  the  transition  point  is  to  be  calculated  according  to  the 

2 

1974  method  of  Granville.  A  second  function  of  body  geometry,  denoted 
by  D,  and  given  by 


(13) 


is  also  calculated  at  each  point  dovmstream  of  the  neutral-stability 
point.  D  will  be  needed  if  the  transition  point  is  to  be  calculated 
according  to  the  1953  method  of  Granville.^ 

CALCULATION  OF  TRANSITION  POINT 

This  section  describes  two  methods  for  calculating  the  point  where 
there  is  a  natural  transition  from  laminar-  to  turbulent-boundary-layer 
flow.  It  Is  to  be  understood  that  natural  transition  may  not  be  reached 
because  forced  transition  takes  place  upstream  of  it,  either  (1)  at  a 
transition  trip  at  a  point  XTRIP  if  ITRIP  is  set  equal  to  1  or  (2)  at  the 

laminar  separation  point  if  ITLS  is  set  equal  to  1. 

Granville  1974  Method 

This  method  computes  Rq  and  TP  at  each  point  downstream  of  the 

neutral-stability  point  and  predicts  transition  at  the  point  where 

first  exceeds  a  particular  polynomial  function  of  TP,  desig- 

nated  by  TF.  The  polynomial  is  selected  by  using  control  variables 

J  and  L  as  input.  J  is  set  equal  to  1  if  and  only  if  it  is  desired  to 

use  a  numerical  approximation  to  the  low-background-turbulence  curve 

2 

shown  as  Figure  16  in  Granville,  together  with  one  of  two  extensions 

of  it  (selected  by  L)  for  large  positive  values  of  TP.  These  are  shown 

in  Figure  3  along  with  the  high-background-turbulence  curve  which  is 
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- LOW  BACKGROUND  TURBULENCt  CURVE.  (FIGURE  16.  PARSONS  ET  AL.' I 

-  -  -  NUMERICAL  APPROXIMATION  TO  LOW  BACKGROUND  TURBULENCE  CURVE 
I  IGURE  16.  PARSONS  E  T  AL  '  U  •  II 


Figure  3  -  Curves  for  Calculating  as  a  Function  of 

and  Numerical  Approximations  to  These  Curves 


selected  If  .1  Is  not  equal  to  1. 
curves  or  TF's  are  (;iven  below. 


The  polynomial  representations  of  these 


J 

L 

TP  Range 

TF 

1 

1  or  /  1 

TP  <  0.075 

671.8  2432  4  *  TP  -  930.4  *  TP^ 

1 

1 

0.75  TP  <  0  165 

2  , 
719.  -  1253.  •  TP  +  5857.9  *  TP  -  7287  7  *  TP^ 

1 

1 

0 165  ^  TP 

639. 

1 

1 

/  1 

0  075  <  TP  <  0.25 

719.  1253.  *  TP  +  5857.9  *  TP^  -  7287  7  *  TP^ 

1 

/  1 

0.25  TP 

551.  +  428.  •  TP 

*  1 

1  or  /  1 

TP  <  0.02 

671  8  -  2432.4  *  TP  -  930.4  *  TP^ 

*  1 

1  or  ^  1 

0.02  <  TP 

681  4  -  1910.38  •  TP  +  1233.6  *  TP^  +  1036  5  '  TP^ 

2 

In  the  concluding  section,  Granville  emphasized  that  the  method  of 
calculating  the  transition  point  is  based  on  empirical  data  for  bodies 
which  do  not  have  significant  adverse  pressure  gradients;  on  a  parallel 
middlebody  TP  is  zero  and  TF  is  either  719.0  or  681.4,  depending  on 
whether  J  is  equal  to  1.  Consequently  the  value  of  Rg  at  transi¬ 
tion,  is  equal  to  plus  a  constant  for  any  transition  point  on  the 

parallel  middlebody.  This  is  too  restrictive  to  be  physically  realistic. 
Nevertheless,  the  method  has  been  used  successfully  with  bodies  which 
have  parallel  middlebodies  in  cases  where  is  large  enough  so  that 
transition  occurs  on  the  forebody.  These  cases  will  be  discussed  later 
in  this  report. 


Granville  1953  Method 

If  it  is  desired  to  use  the  method  of  calculating  the  transition 
point  described  by  Granville  in  1953,^  the  control  variable  153  should 
be  set  equal  to  1.  The  method  predicts  transition  at  the  first  point 
where  the  equivalent  two-dimensional  Rg,  designated  by  superscript  , 

exceeds  its  value  at  the  neutral-stability  point  plus  450  +  400  e^^ 
where  A  is  given  by 


15 


A  - 


1 

5 


4 

9 


(Rg  y/L)^ 

14. 


(Rg<">  y<">/L)^  ' 

R], 


(14) 


The  superscript  n  means  that  quantities  are  to  be  evaluated  at  the 
neutral-stability  point,  ind  D  is  given  by  Equation  (13).  The  equivalent 
two-dimensional  is  given  by 


where  D  is  the  maximum  diameter.  This  method  of  predicting  transi¬ 
tion  can  be  used  on  parallel  mlddlebodies. 

The  form  in  which  this  transition-prediction  method  is  used  is  one 

3  1 

suggested  by  Smith  based  on  a  curve  fit  used  by  White  (see  Equation 

(5-51)  in  White)  to  describe  data  presented  in  Equation  (53)  and  Figure  4 

of  Granville .  ^ 

CALCULATION  OF  TURBULENT  BOUNDARY  LAYER 

The  method  by  which  the  turbulent  boundary  layer  is  calculated  is 
described  in  Granville.  It  is  a  two-parameter  similarity-law  integral 
method  which  does  not  require  the  boundary  layer  to  be  thin  compared  with 
the  body  radius,  and  thus  it  is  valid  near  the  tail  of  an  axisymmetric 
body.  The  equations  needed  for  the  calculations  are  repeated  below.  It 
should  be  emphasized  that  the  method  assumes  that  0  (and  hence  R„)  is 
continuous  through  the  transition  point. 

the  shape  parameter  at  transition,  is  found  from  the  flat- 
plate  equation 


H 


(t) 


1 


0.617  +  0.3863  log  R. 

e  0 


(15) 


13 

White,  F.  M. ,  "Viscous  Fluid  Flow,"  McGraw-Hill  Book  Company, 
New  York  (1974),  pp.  441-444. 
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At  points  dovmstream  from  transition,  o  (defined  by  o  ■  /2/Cj.  where 
C  is  the  local  skin-friction  coefficient)  is  found  from  the  transcendental 


equat ion 


Q_.jA62  (3.889  -Jl)  ^  2.448  log  a  =  2.606  log  R„  -  1.^56 
H  e  e  D 

(H  -  1)0-9392 

-  2.606  log^  1.9392-- 
H 


G,  the  Rotta  shape  parameter,  is  given  by 


G  *  a 


H  -  1 


H,  the  entrainment  shape  parameter  is  given  by 


H  /,  .  1.235  .  33.96  \ 


(e ) 

and  P  is  given  by 


.1.1 


After  the  above  quantities  have  been  found  at  a  given  point,  the 

/s 

reduced  entrainment  factor  E  is  found  from 


3  log^  Rq 


\  ^  H  + 

7  1  +  - 


h  (H  +  1)  ^ 

G  H 

2.606  (H  -  1)‘ 
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E,  the  entrainment  factor,  la  given  by 


E 


(21) 


H,,  the  quadratic  momentum-shape  parameter.  Is  given  by 
<P 


H 


4) 


0.1028  (H  +  3.336) 
H  -  1 


1  0.4746 
o  H  -  1 


1_  7.818 

3  4 

0^  (H  -  1)^ 


(22) 


and 


the  quadratic  displacement-shape  parameter  H 


A 


is  found 


f  rom 


H 


A 


0.4457 
H  -  1 


7.818 
(H  -  1)^ 


(23) 


The  two  boundary-layer  parameters  H  and  Rq  (or  log^  Rg)  play  the 
roles  of  Independent  variables  in  these  calculations.  Derivatives  of 
logg  o,  H,  and  with  respect  to  H  and  log^  Rg  will  be  needed  in  the 
calculations.  It  is  seen  from  Equations  (16),  (18),  (22),  and  (23)  that 
these  are  given  by  the  following  equations: 


i  ^  _  2.8885  1.3469  (H  -  l)o  +  2.606  H(H  -  1.9392) 

o  3H  “  H(H  -  1)  (3.889  -  H)o  +  7.07  H 


(24) 


1  9o  ^ _ 7.527  H _ 

0  9  (log  R„)  (3.889  -  H)o  +  7.07  H 

e  y 


(25) 


9H 

9H 


H(H  -  2)  -  H 
H(H  -  1) 


/  1. 235  ,  93.39  \  P  1  .  ^  1 

H  -  1  V  G  ^2.75  /  L  H(H  -  1)  o  9H  J 


(26) 
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(27) 


3H 


/  1.235  93.39  \  1  da 

H-1  \  G  (.2.75;  0  9  (log^  Rg) 


(})  _  0.1028  H  (2  H  +  0.336  H  -  6.672)  2.606  H’’ 

9H  ■  -  H-1 


0.1821 


+  9  r  »  1  M  i  .  0. 

L o(H  - 1)  J  ( o  is  ^ 


4746  H  (2  H  -  3) 
(H  -  1)^  0 


15.636  (H  -  3)  /  H  \ 

„3  \H  -  1  } 


(28) 


^^4)  r  0.4746  23.45  1  1  9o 

a  (log^  Rq)  “  -  L  (H  -  1)0  -  I)"*  -i  ^ 


(29) 


=  0.4457  (2  H  - 


23.45 
(H  -  1)^ 


1  9o 
3  0  9H 

0 


+ 


U  _  T 

15.636 

0 


5 


(30) 


and 


9 


23.45  .  i  9o  _ 

H-1  ^309  (log^  Rq) 


(31) 


The  order  in  which  these  equations  are  solved  at  each  point  down¬ 
stream  of  transition  is  as  follows.  With  H  and  Rg  known,  Equation  (16) 
is  solved  for  a.  Then,  in  order.  Equation  (17)  is  solved  for  G,  Equa¬ 
tion  (18)  is  solved  for  H,  Equation  (24)  is  solved  for  ^  Equation  (25) 
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Is  solved  for 


Is  solved  for 


i  ao 

0  9  (log^  Rg) 
3H 


Equation  (26)  Is  solved  for 

Equation  (19)  is  solved  for 


Equation  (27) 

Equation  (20) 


is  solved  for  E,  Equation  (21)  is  solved  for  E,  Equation  (22)  is  solved 
for  Equation  (23)  is  solved  for  and  Equations  (28)-(31)  are  solved 
for  the  derivatives  of  and  H^.  Then,  with  superscript  (t  +  1)  denoting 
quantities  at  the  first  point  dovmstream  of  transition,  the  following 
equations  are  used  to  find  Q  and  at  that  point: 


fl(t^l)  /^^l)  9(t^l) 


(t+i) 


± 


(sec  a) 


(t+1) 


,(t+l) 


(32) 


and 


(t+1) 


=  1 


(t+1) 


H 


(t+1)  e 


(t+1) 


(sec  a) 


(t+1) 


(t+1) 


2 


2 


(33) 


At  all  points  downstream  of  transition.  A*  is  found  from 


2  L  L  sec  a  \  L  / 

Li 


(34) 


At  points  downstream  of  the  point  immediately  downstream  of  transi¬ 
tion,  9.  and  are  found  by  integration  of 


/Q ) 

,  y  sec  a  A*/L^  +  2  Q/L^  d  j 

f^L] 

[l^J 

L  2  U/U  ds  ' 

0  °° 

(35) 


20 


and 


21 


■^[h 


1+  ^  i 

rj  sec  a  L 


(H  +  H)' 


+  (H  +  H) 


3  (logg  R0)  3  (iogg  Rq] 


(!)' [<-")  (-t)- 


+  z8 _ an 

"  L  ^  U/  L  L  3  (log  R.)  U/U 

®e  0  “> 


(i)^ 

sec  a  \  L / 


(H  +  H) 


9  doge 


3  (log 


A 

e  “e). 


A(U/U^) 


Equations  (35)-(38)  enable  and  H  to  be  found 


at  the 


(i  -  th  +  1)  point  in  terras  of  functions  evaluated  at  the  i-th  point  The 
values  at  the  (i  -  th  +  1)  point  enable  Equations  (16)-(31)  and  Equa¬ 
tion  (34)  to  be  used  to  evaluate  the  remaining  functions,  so  that  the 
process  can  be  continued  to  the  end  of  the  body. 

CALCULATION  OF  DRAG 

Equation  (7)  in  Granville^  gives  the  formula  for  drag  coefficier.t 
based  on  a  reference  area  A  as 


C  =  4  IT  - j- 

^  A/L^ 
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where  Is  the  momentum  area  of  the  wake  far  downstream.  This  equa¬ 

tion  Is  the  result  of  considering  the  momentum  balance  of  a  flow  through 
a  control  volume  which  contains  the  body  and  has  dimensions  sufficiently 
large  that  the  pressure  at  all  points  on  Its  surface  Is  effectively 
equal  to  the  pressure  In  the  undisturbed  flow. 

If  the  superscript  (e)  denotes  values  at  the  tall  of  the  body, 
can  be  found  from  Equation  (95)  In  Granville^  (with  q  *  7,  as  suggested) 
as 


a 


(D) 


L 


2 


(48) 


where  h  Is  the  ax 1 symmetric  shape  parameter  defined  by 


h  =  A*/S1 


Therefore 


Is  given  by 


h 


(e) 


(49) 


(50) 


COMPARISON  OF  PREDICTED  AND  MEASURED  DRAG 
COEFFICIENTS  AND  TRANSITION  POINTS 

Data  on  bodies  of  revolution  with  and  without  parallel  middlebodies 

are  available  from  recent  experiments.  The  body  without  parallel  middle- 

body  is  represented  by  the  DOLPHIN,  for  which  measurements  were  made 

14  15 

during  drop  tests  in  the  Pacific  Ocean.  ’  Data  on  three  bodies  of 
revolution  with  parallel  middlebodies  are  available  from  recent  DTNSRDC 


^^Carmichael ,  B.  H. ,  "Underwater  Drag  Reduction  through  Choice  of 
Shape,"  American  Institute  of  Aeronautics  and  Astronautics  Paper  66-657 
(1966)  . 

^^Carmlchael ,  B.  H. ,  "Underwater  Drag  Reduction  through  Optimum  Shape," 
in  "Underwater  Missile  Propulsion,"  edited  by  Leonard  Greiner,  Compass 
Publications,  Arlington,  Virginia  (1967),  pp.  147-169. 


23 


towiriR-taiik  tests  of  Models  4620-2,  4620-3,  and  4620-4.  The  speed  range 
of  the  towing-tank  experiments  included  the  range  where  the  length  Froude 
number  is  between  0.5  and  0.6,  in  which  the  theoretical  maximum  of  the 
wave-drag  coecficients  of  submerged  streamlined  bodies  of  revolution 
occur.  No  attempt  has  been  made  to  correct  for  wave  drag,  but  this 
Froude-number  range  is  indicated  on  figures  which  show  the  results  of 
the  towing-tank  experiments.  In  all  of  these  figures,  the  Schoenherr 
friction  line  for  a  flat  plate  with  turbulent  boundary  layer  is  drawn 
in  for  reference  since  the  calculation  method  reproduces  this  line  if 
the  pressure  gradient  is  set  equal  to  zero,  the  body  radius  is  set  equal 
to  a  sufficiently  large  constant,  and  transition  is  forced  at  the  nose. 

In  this  section,  is  the  drag  coefficient  based  on  wetted  area, 
x^  is  the  axial  distance  from  the  nose  to  the  transition  point,  L  is  the 
model  length,  and  and  F^  are  the  Reynolds  number  and  Froude  number 
based  on  L  and  model  speed. 

DOLPHIN 

The  DOLPHIN  profile  is  that  of  an  NACA  66-030  airfoil  section, 

faired  into  a  boom  at  the  stern  on  which  stabilizing  fins  are  mounted. 

Drop  tests  were  performed  in  the  Pacific  Ocean  and  speed  was  computed 

from  a  time  history  of  the  dynamic  pressure  at  the  model  nose.  The 

reporteii  drag  is  for  the  "airfoil"  portion  of  the  model  only;  the  boom 

14,15 

and  fin  drag  were  estimated  by  a  method  described  by  Carmichael 
and  subtracted  out.  Consequently  the  drag  coefficients  reported  here 
are  for  a  model  truncated  at  the  stern  of  the  airfoil  profile.  The 
actual  DOLPHIN  and  the  truncated  approximation  are  shown  in  Figure  4. 

The  offsets  of  the  approximation  are  listed  in  Table  1. 

Figure  5  shows  the  predicted  drag  coefficients  and  transition 
points,  the  "measured"  drag  coefficients,  and  transition  points  deduced 
from  the  measured  drag  coefficients  according  to  the  method  of  Young. 
Considerable  scatter  among  the  results  is  apparent;  this  might  be 


Young,  A.  D.,  "The  Calculation  of  the  Total  and  .Skin  Friction  Dr.igs 
of  Bodies  of  Revolution  at  Zero  Incidence,"  British  Aeronautical  Researcli 
Committee  Reports  and  Memoranda  1874  (Apr  1939). 


I 
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Figure  5  -  Predicted  and  Measured  Cp  and  for  DOLPHIN 


TABLE  1  -  OFFETS  AND  DIMENSIONS  OF  APPROXIMATED  DOLPHIN 


X 

L 

Y 

L 

X 

L 

Y 

L 

0.00 

0.00000 

0.35 

0.14565 

0.01 

002768 

0.40 

0.14894 

0.02 

0.03920 

0.45 

0.15000 

0.03 

0.04808 

0.50 

0.14874 

0.04 

0.05556 

0.55 

0.14434 

0.05 

0.06215 

0.60 

0.13635 

0.06 

0.06810 

0.65 

0.12478 

0.07 

0.07354 

0.70 

0.10992 

0.08 

007859 

0.75 

0.09228 

0.09 

0.08329 

0.80 

0.07262 

0.10 

0.08770 

0.85 

0.05187 

0.15 

0.10635 

0.90 

0.03133 

0.20 

0.12072 

0.95 

0.01287 

0.25 

0.13178 

1.00 

0.00000 

0.30 

0.14000 

Total  Length,  ft  6.300 

Wetted  Area,  ft^  18.150 

Volume,  ft^  5.600 


attributed  to  variations  in  the  ocean  ambient-turbulence  levels  between 
tests  since  such  variations  are  believed  to  affect  the  mechanism  of 
boundary-layer  transition.  Nonetheless,  the  lines  representing  drag 
coefficients  and  transition  points  based  on  the  1953  prediction  of 
transition  are  seen  to  lie  closer  to  the  measured  drag  coefficients  and 
deduced  transition  points  than  do  the  lines  of  the  1974  transition 
prediction. 


MODELS  4620-2,  4620-3,  AND  4620-4 

These  models  have  streamlined  forebodies  and  afterbodies  separated 
by  parallel  middlebodies .  Note  from  the  model  offsets  (Table  2)  that 
Model  4620-2  has  the  bluntest  nose  of  the  three  and  that  Model  4620-4 
has  the  finest  nose.  The  models  were  towed  by  vertical  struts  extending 
down  from  the  towing  carriage  into  the  tops  of  the  models  near  amidships. 
Each  model  was  tested  in  a  bare-hull  condition  and  with  an  0.024-inch 


TABLE  2 


OFFSETS  AND  DIMENSIONS  OF  MODELS  4620-2,  4620-3,  AND  4620-4 


X 

Y 

Y 

Y 

— 

— .  4620  2 

—  .  4620  3 

— ,  4620  4 

L 

L 

L 

L 

0.00 

0  00000 

0.00000 

0  00000 

0.01 

0  01954 

0.01431 

001096 

0.02 

0.02753 

0.02023 

0.01549 

0.03 

0.03339 

0.02472 

0.01897 

0.04 

0  03789 

002845 

0.02188 

0  05 

0  04127 

0.03163 

002443 

0.06 

0.04363 

0.03437 

0.02670 

0  07 

0  04507 

0.03674 

002875 

0.08 

0  04573 

0.03876 

0.03062 

0.09 

0.04045 

003233 

0.10 

0.04183 

0.03388 

0.15 

0.04472 

0.03967 

0.20 

P.M 

.B 

1 

0.04254 

0  25 

from 

0.04322 

0.30 

xa 

0.088 

A 

to 

P.M. 8. 

1 

0.35 

x'V 

0.60 

ftom 

P.M.B. 

0.40 

with 

xA 

0.16 

from 

0.45 

0.50 

vA'  0.04587 

TO 

xA  =  0.26 

xA  0.61 

with 

to  ^ 

x,/l'  -  0.62 

0.55 

' 

yA  =  0 

04475 

with 

yA  =  0.04323 

\ 

0.60 

0.04587 

0.65 

0.04543 

0.04449 

0.04312 

0.70 

0.04345 

0.04280 

0.04182 

0.75 

0.03986 

0  03942 

0.03874 

0.80 

0.03502 

0.03471 

0.03421 

0.85 

0.02915 

0.02892 

0.02856 

0.90 

0.022C't 

0.02186 

0.02164 

0.95 

0.01273 

0.01268 

0.01261 

1.00 

0.00000 

0.00000 

0.00000 

Models 

4620  2 

4620  3 

4620  4 

Tnt.'il  Lenqth.  It 

22  300 

22.859 

23  664 

2 

Wetted  Aren,  ft 

123,871 

124  764 

126  345 

3 

Volume,  ft 

58.102 

58  102 

58.102 
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wire  ring  mounted  at  5  percent  of  the  axial  distance  aft  of  the  nose  to 
stimulate  boundary-layer  transition  at  that  station.  In  the  presentation 
of  the  drag  coefficients  measured  In  these  experiments,  no  attempt  has 
been  made  to  estimate  the  drag  on  the  wire  ring  and  subtract  it  out. 

The  experiments  conducted  with  these  models  are  described  in  detail 
by  McCarthy,  Power,  and  Huang. 

The  drag-coefficient  results  of  Model  4620-2  are  shown  in  Figure  6. 
This  model  is  so  blunt  that  the  laminar-separation  point  was  only  8.1 
percent  of  the  axial  length  aft  of  the  nose;  both  methods  predict  that 
transition  in  the  bare-hull  experiments  will  take  place  there  throughout 
the  speed  range  used.  It  is  seen  that  below  the  Reynolds-number  range 
where  there  is  appreciable  wavemaking  drag,  the  drag  prediction  is  fairly 
good;  the  overprediction  is  generally  no  more  than  3  percent  for  both 
artificially  tripped  transition  and  transition  which  takes  place  at  the 
laminar-separation  point.  However,  the  overprediction  in  the  former 
case  is  actually  larger  because  trip-wire  drag  is  included  in  the 
measurements . 

The  drag  coefficients  and  transition  points  for  Model  4620-3  are 
shown  in  Figure  7.  The  1953  transition-point  prediction  method  is  accu¬ 
rate  except  at  the  lower  Reynolds  numbers,  where  both  methods  predict 
transition  much  further  back  on  the  parallel  middlebody  than  the  point 
at  which  it  was  measured.*  At  Reynolds  numbers  where  transition  occurred 
on  the  forebody,  the  1953  method  predicted  its  location  to  within  about 
6  percent;  the  1974  method  predicted  it  slightly  further  aft.  In  the 
Reynolds-number  range  between  about  12  and  22  million  (for  which  the 
1953  method  predicted  the  transition  point  to  this  accuracy  and  wave¬ 
making  drag  was  negligible)  the  drag  coefficient  predicted  by  using  the 
1953  method  is  quite  accurate;  it  overpredicts  the  drag  by  less  than  2 
percent.  However,  drag  overprediction  is  larger  with  tripped  transi¬ 
tion;  5  percent  is  representative.  Again,  this  is  actually  larger 
because  the  measured  drag  includes  the  drag  of  the  trip  wire. 

^^McCarthy,  J.  H.  et  al.,  "The  Roles  of  Transition,  Laminar  Separation, 
and  Turbulence  Stimulation  in  the  Analysis  of  Axisymmetrlc  Body  Drag," 
DTNSRDC  Report  4728  (in  review). 

7 

*The  Granville  1974  correlation  does  not  include  bodies  with  parallel 
middlebodies  and  hence  is  not  strictly  applicable  here. 


Figure  8  shows  the  results  for  Model  4620-4.  The  1953  method  pre¬ 
dicted  the  transition  point  to  within  7  percent  in  the  Reynolds-number 
range  between  about  10  and  17  million;  at  higher  Reynolds  numbers,  the 
1974  method  was  tlie  more  accurate,  predicting  it  to  within  10  penent; 
below  this  range,  botli  methods  predicted  transition  much  further  aft  of 
where  it  was  measured.  In  the  Reynolds-number  range  between  about  10  and 
20  million,  where  the  1953  method  is  fairly  accurate,  the  drag  coefficient 
is  predicted  to  within  3  percent.  The  1974  method  predicted  the  transi¬ 
tion  point  more  accurately  at  Reynolds  numbers  above  this  range;  it  also 
predicted  the  drag  ciicfficient  to  within  3  percent  at  Reynolds  numbers 
between  20  and  25  million,  above  which  wavemaking  drag  becomes  appre¬ 
ciable.  The  drag  coefficient  with  tripped  transition  is  overpredicted; 
once  more,  the  discrepancy  is  actually  greater  than  it  appears  because 
trip-wire  drag  is  included  in  measured  drag. 

DiSCUSSION  AND  CONCLUSIONS 

The  most  important  conclusion  to  be  drawn  from  the  foregoing  is 

4 

that  the  updated  turbulent-boundary-layer  theory  of  Granville  is 
capable  of  giving  reasonable  predictions  of  the  viscous  drag  on  stream¬ 
lined  bodies  of  revolution  as  long  as  the  transition  point  is  located 
accurately.  The  digital-computer  program  makes  it  possible  to  use  the 
theory  on  a  routine  basis;  no  numerical  difficulties  have  arisen.  In 
cases  where  the  transition  point  is  fixed  by  tripping  the  laminar 
boundary  layer  or  where  the  transition  point  is  predicted  accurately, 
comparison  with  measured  drag  shows  the  theory  and  computer  program 
consistently  overpredict  the  drag  coefficient  by  as  much  as  5  percent. 

This  may  be  due  to  the  values  of  the  boundary-layer  constants  used 
which  have  conriderable  experimental  variation. 

Unfortunately,  two  factors  which  have  been  neglected  would  each 
increase  the  overprediction  had  they  been  considered.  First,  the  addi¬ 
tional  drag  due  to  the  presence  of  displacement  thickness  is  not  present 
in  the  predictions  since,  in  the  cases  presented  here,  the  drag  calcula¬ 
tions  were  not  repeated  after  the  appropriate  modification  to  the 
hydrodynamic  source  strength.  Second,  in  cases  where  the  trip  wire  was 


32 


mounted,  the  measured  drag  includes  the  drag  on  the  wire  and  thus  is 
artificially  increased.  Furthermore,  the  overprediction  was  present 
(and  often  at  its  worst)  when  transition  was  fixed  at  only  5  percent  aft 
of  the  nose.  Thus  it  cannot  be  attributed  to  an  inaccuracy  in  the 
laminar-boundary-layer  calculations.  As  was  mentioned,  the  program  re¬ 
produces  the  Schoenherr  frictional  line  for  a  flat  plate  with  turbulent 
boundary  layer  under  the  appropriate  circumstances.  Therefore,  the 
overpredJ ctlon  should  be  attributed  to  the  means  by  which  the  presence 
of  a  pressure  gradient  is  Included  in  the  turbulent-boundary-layer 
calculations,  either  on  the  body  itself  or  in  the  turbulent  wake  since 
the  drag  calculation  is  based  on  a  particular  model  of  the  wake. 

A  better  means  of  estimating  the  location  of  the  transition  point 
is  needed  since,  at  Reynolds  numbers  of  10  million  or  less,  an  inaccurate 
prediction  of  the  transition  point  results  in  a  serious  error  in  the 
drag  prediction;  see  Figures  7  and  8.  Comparison  between  the  two  transition- 
point  prediction  methods  is  inconclusive.  Figures  7  and  8  show  tnat 
neither  is  accurate  if  transition  occurs  on  the  parallel  middlebody, 
although  the  1974  method  is  not  intended  for  this  case.  Aside  from 
this.  Figures  5,  7,  and  8  show  that  the  1953  method  was  more  accurate 
for  DOLPHIN,  Model  4620-3,  and  the  lower  part  of  the  Reynolds-number 
range  for  Model  4620-4,  whereas  the  1974  method  was  more  accurate  in  the 
upper  part  of  the  Reynolds-number  range  for  Model  4620-4. 

Comparison  of  the  location  of  the  Schoenherr  friction  line  for  a 
flat-plate  turbulent  boundary  layer  relative  to  the  measured  drag 
coefficients  in  Figure  5  with  its  same  relative  location  in  Figures  6, 

7,  and  8  shows  that  the  designers  of  DOLPHIN  were  successful  in  selecting 
a  body  of  revolution  which  maintained  an  appreciable  extent  of  laminar 
boundary  layer  throughout  a  practical  range  of  Reynolds  numbers. 
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DESCRIPTION,  GLOSSARY,  LISTING,  AND  SAMPLE 
RUN  OF  DP INI 


DESCRIPTION  OF  DP  INI 

18 

This  program  uses  the  method  of  Granville  to  compute  the  offsets 

and  slopes  of  a  streamlined  body  of  revolution  formed  by  a  forebody  with 

rounded  nose,  no  parallel  middlebody,  and  an  afterbody  with  pointed 

tail.  The  axis  of  symmetry  is  the  x-axis,  the  nose  is  at  the  origin, 

and  the  tail  is  at  x  *  1.0.  On  the  forebody,  y  is  equal  to  the  square 

root  of  a  quartic  in  x;  on  the  afterbody,  it  is  equal  to  the  square  root 

of  a  quintic.  Coefficients  in  the  quartic  and  quintic  are  specified  by 

geometric  parameters  which  are  Input  variables  together  with  the  require- 

2  2 

ment  that  the  slope  dy/dx  be  zero  and  that  the  curvature  d  y/dx  be 
continuous  at  the  £orebody-af terbody  juncture,  which  is  the  location  of 
maximum  diameter. 

As  calculated  by  this  program,  the  bodies  are  represented  by  series 
of  frustums  of  cones  with  generating  axis  along  the  axis  of  revolution 
since  the  body  shapes  are  formed  by  straight-line  segments  between 
points  defined  by  pairs  of  x-  and  y-  coordinates.  One  feature  of  the 
program  is  that  the  segment  lengths  of  step  sizes  over  the  foremost 
portion  of  the  body  are  one-tenth  as  large  as  those  over  the  rest  of  it. 
(A  different  ratio  of  step  sizes  could  be  attained  by  a  minor  modifica¬ 
tion  to  the  appropriate  DPIN  and  to  DPOUT.)  The  change  from  the  smaller 
to  the  larger  step  size  is  made  over  a  total  of  three  steps  in  such  a 
way  th£)t  no  step  is  more  than  twice  as  large  as  the  one  preceding  it. 
Experience  in  operation  of  the  axisymmetric  potential-flow  program  has 
shown  the  desirability  of  this  relatively  gradual  change. 

Input  variables  to  the  program  consist  of  the  integers  NPTS,  NCPTS, 
and  INFL  and  the  floating-point  numbers  M,  KL,  RL,  SL,  and  ELD.  NPTS 
and  NCPTS  are  the  total  number  of  points  on  the  body  and  the  number  of 
closely  spaced  points.  INFL  is  the  maximum  number  of  inflection  points 


18 

Granville,  P.  S.,  "Geometrical  Characteristics  of  Streamlined  Shapes," 
NSRDC  Report  2962  (Mar  1969);  also  Journal  of  Ship  Research,  Vol.  13, 

No.  4,  pp.  299-313  (Dec  1969). 
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permitted  on  the  afterbody;  it  must  be  either  0  or  1.  M  is  the  ratio  of 
the  forebody  length  to  total  length,  KL  is  the  curvature  at  maximum 

diameter,  RL  is  the  radius  of  curvature  at  the  nose,  SL  is  the  slope  at 

the  tail,  and  ELD  is  the  ratio  of  total  length  to  maximum  diameter.  KL, 

RL,  and  SL  are  expressed  in  coordinates  scaled  so  that  the  total  length 

of  the  body  is  1. 

18 

In  order  to  use  the  Granville  polynomials,  they  are  transformed 
to  a  stretched  coordinate  in  which  the  forebody  or  afterbody  length  and 
the  maximum  radius  each  equals  1.  Thus  the  coef f iclerits  in  the  quartic 
for  the  forebody  are  R,  the  stretched  nose  radius,  and  KF,  the  stretched 
curvature  at  maximum  diameter.  The  coefficients  in  the  quintic  for  the 
afterbody  are  KA,  the  (differently)  stretched  curvature  at  maximum  diam¬ 
eter,  and  S2,  the  square  of  SP,  the  stretched  slope  at  the  tail.  KF, 

K,^,  and  SP  also  have  their  signs  reversed  so  as  to  be  ordinarily  posi- 

1 8 

tive.  The  Granville  method  of  avoiding  undesirable  bulges,  zeros,  and 
inflection  points  on  the  forebody  requires  R  and  KF  to  lie  in  an  admiss¬ 
ible  region  when  graphed  (see  Figures  2-4  of  Reference  18,  in  which  they 

are  denoted  by  r  and  k^).  In  the  same  way,  undesirable  afterbodies  are 
excluded  by  requiring  that  S2  and  KA  lie  in  one  of  two  admissible  regions 

when  graphed  (see  Figures  6  and  7  of  Reference  18,  in  which  they  are 
2 

denoted  by  s  and  k^^),  according  to  whether  zero  or  one  inflection 
points  are  to  be  allowed  on  the  afterbody. 

Instead  of  introducing  the  complications  associated  with  numerical 
representation  of  the  curves  which  delineate  the  admissible  regions  on 
the  graphs,  DPINl  uses  a  different  method  for  checking  the  occurrence  of 
bulges  and  points  of  inflection:  as  the  slope  at  each  point  is  computed, 
it  is  compared  with  zero  and  the  slope  at  the  previous  point.  If  at 
points  on  the  forebody  the  slope  is  negative  or  greater  than  the  previous 
slope,  the  program  prints  out  a  message  that  the  forebody  parameters  are 
inadmissible  and  stops.  At  points  on  the  afterbody,  the  program  prints 
out  the  corresponding  message  and  stops  if  the  slope  is  positive  or,  if 
no  inflection  points  are  permitted  (INFL  =0),  the  slope  is  greater  than 
the  previous  slope.  If  one  inflection  point  is  permitted  (INFL  =  1), 
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the  slope  Is  allowed  to  exceed  the  previous  slope  over  a  range  of  points, 
past  which  It  may  become  less  than  the  previous  slope  over  one  range  but 
must  not  exceed  the  previous  slope  over  a  second  range. 

The  program  user  may  wish  to  modify  DPINl  so  that  KF,  KA,  R,  and  S2 

(or  Its  square  root)  are  used  as  input  rather  than  the  unstretched  KL, 

RL,  and  SL.  Furthermore,  he  may  wish  to  replace  either  KF  or  R  with 

CPF,  the  forebody  prismatic  coefficient,  or  to  replace  either  KA  or  S2 

with  CPA,  the  afterbody  prismatic  coefficient,  since  Equations  (265)  and 
18 

(266)  in  Granville  present  linear  relations  between  these  respective 
triads  of  geometric  parameters. 

After  DPINl  has  read  its  seven  input  variables,  it  writes  them  out. 
Next  it  computes  the  stretched  variables  and  the  forebody,  afterbody, 
and  overall  prismatic  coefficients  and  writes  them  out.  Then  it  finds 
the  appropriate  step  size  DX  and  proceeds  to  calculate  the  offsets  Y(I) 
at  each  point  along  the  forebody  and  afterbody  by  taking  the  square  root 
of  the  appropriate  polynomial.  The  slope  DYDX(I)  and  increment  in  arc 
length  DS  are  calculated  by  using  the  derivative  of  the  polynomial.  The 
wetted  area  and  a  geometric  function  TP  (I),  which  is  used  for  predicting 
boundary-layer  transition,  are  also  computed  at  each  point.  Then  I,  X, 

Y,  DYDX,  S,  and  TP  are  written  out  for  each  point.  Finally  the  total 
volume,  forebody  and  afterbody  lengths,  and  the  forebody,  afterbody,  and 
overall  wetted  areas  and  arc  lengths  are  written  out.  X(I),  Y(I), 
DYDX(I),  and  the  volume,  wetted  area,  length-to-diameter  ratio,  and 
numbers  of  points  and  of  close  points  are  also  written  on  tape  for  use 
by  subsequent  programs. 

DPINl  -  GLOSSARY 

The  following  glossary  of  variables  used  in  DPINl  is  arranged 
alphabetically  by  FORTRAN  variable  name. 


39 


FORTRAN 

Variable 

Name 

Variable 

Definition 

AA 

a("Vl2 

Wetted  area  of  afterbody  divided  by  length 
squared 

AF 

Wetted  area  of  forebody  divided  by  length 
squared 

AT 

Total  wetted  area  divided  by  length  squared 

CP 

S 

Prismatic  coefficient  of  complete  body 

CPA 

c 

P 

Prismatic  coefficient  of  afterbody 

CPF 

c 

p 

Prismatic  coefficient  of  forebody 

DS 

ds/L 

Increment  in  arc  length  s  divided  by 
total  length 

DX 

dx/L 

Increment  in  axial  length  x  divided  by 
total  length 

DYDX 

dx 

Slope,  dimensioned  to  be  a  function  of  I 

ELD 

L/MAX.  DIAMETER 

Ratio  of  total  length  to  maximum  diameter 

I 

i 

Integer  which  increases  from  1  at  the  nose 
to  NCPTS  at  the  point  where  step  size 
begins  increasing,  thence  to  NPTS  at  the 
tail 

INFL 

Control  variable  which  specifies  the  maxi¬ 
mum  number  of  inflection  points  allowed 
on  the  afterbody;  always  input  as  0  or  1 

KA 

k 

Curvature  at  forebody-afterbody  juncture, 
normalized  for  unit  afterbody  length  and 

unit  maximum  radius,  with  sign  reversed 
so  as  to  be  ordinarily  positive 
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FORTRAN 

Variable 

Name 

Variable 

Definition 

KF 

Curvature  at  forebody-afterbody  juncture, 
normalized  for  unit  forebody  length  and 
unit  maximum  radius,  with  sign  reversed 
so  as  to  be  ordinarily  positive 

KL 

L  ^ 
dx^ 

(f orebody- 
af  terbody 
j  uncture) 

Total  length  times  curvature  at  forebody- 
afterbody  juncture 

M 

m 

Ratio  of  forebody  length  to  total  length 

Ml 

1  -  m 

Ratio  of  afterbody  length  to  total  length 

NCPTS 

Number  of  the  point  downstream  of  which 
the  point  spacing  begins  to  increase 

NCI 

n(^)  4 

•  1 

NC2 

n(^)  4 

■  2 

NC3 

n(^)  4 

■  3 

NC4 

n<^>  4 

■  4 

NPTS 

n 

Total  number  of  points 

PI 

TI 

Ratio  of  the  circumference  of  a  circle 
to  its  diameter 

R 

r 

Radius  of  curvature  at  nose,  normalized 
for  unit  forebody  length  and  unit  maxi¬ 
mum  radius 

RL 

1  d^x 

L  ^  2 

dy 

(nose) 

Radius  of  curvature  at  nose,  divided  by 
total  length 

S 

s/L 

Arc  length  at  arbitrary  x,  divided  by 

total  length;  dimensioned  to  be  a 
function  of  1 
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FORTRAN 

Variable  Variable 
Name 


SA 


SF 


Def initlon 


Arc  length  of  afterbody  divided  by  total 
length 

Arc  length  of  forebody  divided  by  total 
length 


SL 


ST 

S2 

TP 

VL3 

X 

XA 


XAl 


XF 


dx 


(tail) 


s 


(P) 


Slope  at  tail 


Slope  at  tail,  normalized  by  unit  after¬ 
body  length  and  unit  maximum  radius,  with 
sign  reversed  so  as  to  be  ordinarily 
positive 

Total  arc  length  divided  by  total  length 


SP  squared 


1  dy 
ELD  y(x)  dx 


Function  used  in  DPOUT  to  predict  transi¬ 
tion;  dimensioned  to  be  a  function  ot  I 


Volume  of  body  divided  by  cube  of  total 
length 


x/L 


Axial  coordinate  divided  by  body  length; 
dimensioned  to  be  a  function  of  I 


1  -  x/L 
1  -  m 


m  -  x/L 
1  -  m 


x/Lm 


Reversed  x,  normalized  to  increase  from 
XA  =  0  at  the  tail  to  XA  =  1  at  the 
forebody-afterbody  juncture 


Reversed  x,  normalized  to  increase  from 
XAl  =  -1  at  the  tail  to  XAl  =  0  at  the 
forebody-afterbody  juncture 

X  normalized  to  increase  from  XF  =  0  at 
the  nose  to  XF  =  1  at  the  forebody- 
afterbody  juncture 
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FORTRAN 

Variable 

Variable 

Definition 

Name 

XFl 

x/L  -  m 
m 

X  normalized  to  increase  from  XFl  -  -1 
at  the  nose  to  XFl  ■  0  at  the  forebody- 
afterbody  Juncture 

XX 

x/L 

Axial  length,  identical  to  x  except  that 
it  is  not  dimensioned  as  a  function  of  I 

Y 

y/L 

Body  radius  divided  by  total  length; 
dimensioned  to  be  a  function  of  I 
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PROGRflMOPINl(lNPUT,CUrPUT,TaPr60,TAPE61,TAPE.5=INPUT,TAPE6  =  OUTPUT,Tft 

UAPE86) 
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CPF  =  .1*R  -KF/30.  » 

CPA  =  (S2-<A) /60. ^.5 
CP  =  M*CPF»M1*CPA 
WRITE  (6,23)  COF,CPA 


23  POPHflT  (/IXt*  CPF  =*,Fll.a,3X,  *CFA  =  •,F11.8,3X,  *CP  =  *f  Fll,8) 
WRITE  (6»13) 

19  PORPAT  (//4X,  *1*,  9X,*X»‘,  1»X,  *V*,  14X,  •OVt)X*,  IIX,  *S*,  14X, 
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26  FORKAK/,  IX,  *F0RE8O0Y  ARC  LENGTH/L  =  *,F10.6,  3X,  •AFTERBODY 
fARC  LEHGTH/L  =  •,F1G.6,3X,  ♦  TOTAL  ARC  LEGNTH/L  =  •,  F10.6) 
HRITE(88)  VL3,AT,  FLO,  NFTS,  NCPIS 


G0T07 
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FORMAK//*  F0RE900Y  PARAMETERS  ARE  INADMISSIBLE*) 
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This  program  uses  the  polynomials  of  Landweber  and  Gertler  to 
compute  the  offsets  and  slopes  of  a  David  Taylor  Model  Basin  Series  38 
body.  This  is  a  streamlined  body  of  revolution  in  which  Y(I)  is  the 
square  root  of  a  septic  in  X(I),  where  the  axis  of  symmetry  is  the 
x-axis,  the  nose  is  at  the  origin,  and  the  tall  is  at  x  ■  1.0.  Coeffi¬ 
cients  in  the  septic  are  specified  by  geometric  parameters.  The  user  is 
warned  that,  in  contrast  to  DPINl  and  DPIN3,  no  checks  are  provided  in 
DPIN2  to  ensure  that  the  geometric  parameters  specified  are  within 
admissible  ranges.  Instead,  it  is  assumed  that  the  user  either  is 
describing  an  existing  Series  38  body  or  has  Reference  19  available  so 
that  the  parameters  have  been  selected  for  him  or  he  can  readily  check 
their  admissibility. 

As  calculated  by  this  program,  the  bodies  are  represented  by  series 
of  frustums  of  cones  with  generating  axis  along  the  axis  of  revolution 
since  the  body  shapes  are  formed  by  straight-line  segments  between 
points  defined  by  pairs  of  x-  and  y-  coordinates.  One  feature  of  the 
program  is  that  the  segment  lengths  or  step  sizes  over  the  foremost 
portion  of  the  body  are  one-tenth  as  large  as  those  over  the  rest  of  it. 
(A  different  ratio  of  step  sizes  could  be  attained  by  a  minor  modifica¬ 
tion  to  the  appropriate  DPIN  and  to  DPOUT.)  The  change  from  the  smaller 
to  the  larger  step  size  is  made  ov^r  a  total  of  three  steps  in  such  a 
way  that  no  step  is  more  than  twice  as  large  as  the  one  preceding  it. 
Experience  in  operation  of  the  axisymmetrlc  potential-flow  program  has 
shown  the  desirability  of  this  relatively  gradual  change. 

Input  variables  to  the  program  consist  of  the  integers  NPTS  and 
NCPTS  and  the  floating-point  numbers  M,  RO,  Rl,  CP,  and  ELD.  NPTS  and 
NCPTS  are  the  total  number  of  points  on  the  body  and  the  number  of 


19 

Landweber,  L.  and  M.  Gertler,  "Mathematical  Formulation  of  Bodies 
of  Revolution,"  David  Taylor  Model  Basin  Report  719  (Sep  1930). 
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closely  spaced  points.  M  Is  the  ratio  of  the  axial  distance  from  nose 
to  maximum  diameter  to  the  total  length.  RO  Is  the  nose  radius  of 
curvature  times  the  total  length  divided  by  the  square  of  the  maximum 
diameter.  R1  Is  the  tail  radius  of  curvature  times  the  total  length 
divided  by  the  square  of  the  maximum  diameter.  CP  is  the  prismatic 
coefficient,  and  ELD  Is  the  ratio  of  total  length  to  maximum  diameter. 

After  DPIN2  has  read  its  seven  input  variables,  it  writes  them  out. 
Then  it  evaluates  the  septic  and  Its  derivative  at  each  point  and  finds 
the  offsets  Y(I)  and  slope  DYDX(I)  in  terms  of  these.  The  increment  in 
arc  length  DS,  the  wetted  area,  and  the  geometric  function  TP(I),  which 
is  used  for  predicting  boundary-layer  transition,  are  also  computed  at 
each  point.  Then  I,  X,  Y,  DYDX,  S,  and  TP  are  written  out  for  each 
point.  Finally  the  total  volume,  wetted  area,  and  arc  length  are  com¬ 
puted  and  written  out.  X(I),  Y(I),  DYDX(I),  and  the  volume,  wetted 
area,  length-to-diameter  ratio,  and  the  numbers  of  total  points  and 
close  points  are  also  written  on  tape  for  use  by  subsequent  programs. 


DPIN2  -  GLOSSARY 

The  following  glossary  of  variables  used  in  DPIN2  is  arranged 
alphabetically  by  FORTRAN  variable  name. 


FORTRAN 

Variable  Variable 
Name 


AT 


AO  a 

o 

A1 

BO  3 

o 

B1 


CP 


c 

p 


Definition 

Total  wetted  area  divided  by  square  of  total 
length 

Coefficient  used  in  calculation  of  ROX 
Coefficient  used  in  calculation  of  ROX 
Coefficient  used  in  calculation  of  RIX 
Coefficient  used  in  calculation  of  RIX 

Prismatic  coefficient 
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FORTRAN 

Variable 

Variable 

Definition 

Name 

DP 

dx 

Derivative  of  PX  with  respect  to  X 

DQ 

dx 

Derivative  of  QX  with  respect  to  X 

dR 

DRO 

o 

dx 

Derivative  of  ROX  with  respect  to  X 

dR 

DRl 

X 

dx 

Derivative  of  RIX  with  respect  to  X 

DS 

ds/L 

Increment  in  arc  length  s  divided  by  total 
length 

DX 

dx/L 

Increment  in  axial  length  x  divided  by  total 
length 

DYDX 

dZ 

dx 

Slope;  dimensioned  to  be  a  function  of  I 

DY2DX 

1  d^ 

L  dx 

2 

Derivative  of  y  with  respect  to  x,  divided 

by  total  length 

DO 

6 

o 

Coefficient  used  in  calculation  of  QX 

D1 

^1 

Coefficient  used  in  calculation  of  QX 

D2 

^2 

Coefficient  used  in  calculation  of  QX 

ELD 

L/MAX.  DIAMETER 

Ratio  of  total  length  to  maximum  diameter 

G 

g 

Coefficient  used  in  calculation  of  PX 

I 

i 

Integer  vhich  increases  from  1  at  the  nose  to 
NCPTS  at  the  point  where  step  size  begins 
increasing,  thence  to  NPTS  at  the  tail 

M 

m 

Ratio  of  forebody  length  to  total  length 

Ml 

1.  -  m 

Ratio  of  afterbody  length  to  total  length 
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FORTRAN 

Variable 

Name 

Variable 

Def inition 

NCPTS 

Number  of  the  point  downstream  of  which 
the  point  spacing  begins  to  Increase 

NCI 

n(^)  + 

1 

NC2 

n(^)  ^ 

2 

NC3 

n(^)  . 

3 

NC4 

4 

NPTS 

n 

Total  number  of  points 

PI 

71 

Ratio  of  the  circumference  of  a  circle  to  its 
diameter 

PX 

P(x) 

Polynomial  in  X,  XI,  and  XM 

QX 

Q(x) 

Polynomial  in  X  and  XI 

RO 

r 

0 

Dimensional  radius  of  curvature  at  nose, 
multiplied  by  body  length  and  divided  by  the 
square  of  the  maximum  diameter 

ROX 

Ro(x) 

Polynomial  in  X,  XI,  and  XM 

R1 

>^1 

Dimensional  radius  of  curvature  of  tail 
multiplied  by  the  body  length  and  divided  by 
the  square  of  the  maximum  diameter 

RIX 

Rj^(x) 

Polynomial  in  X,  XI,  and  XM 

S 

s/L 

Arc  length  at  arbitrary  X  divided  by  total 

length;  dimensioned  to  be  a  function  of  I 


ST 

s(^)/L 

Total  arc  length  divided  by  total 

length 

TP 

L 

ELD  y(x) 

dx 

Function  used  in  DPOUT  to  predict 
dimensioned  to  be  a  function  of  I 

transition; 

VL3 

V/L^ 

Volume  of  body  divided  by  cube  of 

total  length 
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FORTRAN 

Variable 

Name 

Variable 

Definition 

X 

x/L 

Axial  coordinate  divided  by  total  length; 
dimensioned  to  be  a  function  of  I 

XM 

X 

--m 

X  stretched  to  Increase  from  -m  at  the  nose  to 
1  -  m  at  the  tall 

XX 

x/L 

Axial  coordinate,  Identical  to  X  except  that 

It  Is  not  dimensioned  as  a  function  of  I 

XI 

i-^ 

X  normalized  to  Increase  from  XI  =  -1  at  the 
nose  to  XI  *  0  at  the  tall 

Y 

y/L 

Body  radius  divided  by  total  length;  dimen¬ 
sioned  to  be  a  function  of  I 

Y2 

y  /L 

Y  squared 
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DO  *  ,75*(2.-8.*M>7,«M*M)/(M*^*M1*M1*M1*(2,-7.*M*7.*M*M) ) 

D1  »  -(1.-2.*M-7,*M*M*10.5*M*M*M)/(M*M*M*M1*M1*M1*(2.-7.*M 
♦♦7.*M*M) ) 

D?  »  1,75*{1.-5.*M4S.*M*M)/(M*M*M*M1*M1*M1*(2.-7.*M^7.*M*M) ) 

DO  2  I=2»NPTS 
IFd.EQ.NCn  OX  «  .002 
ird  .EG.  N02)  DX  =  .003 
IFd  .EO.  NC3)  OX  =  .005 
IFd  .EQ.  NC4)  OX  =  .010 
Xd)=X(I-l)^OX 
XX  =  X(I) 

XI  =  XX-1. 

XM  *  XX-M 

ROX  =  XX*X1*X1*XM*XM*(A0^A1*XX) 

RIX  =  XX*XX*X1*XM*XM* (flO*Bl*XX) 

PX  =  XX*XX*X1*X1*XM*XM*G 

OX  =  XX*XX*X1^X1*(00>D1*XX4D2*XX*XX) 

Y2  »  2.*R0*W0X^2.*Rl*RlX^CP*PX»ex 
Yd)  =  SQRT(Y2)/ELn 

DRO  =  R0X«(l./XX^2./Xl^2./XM^Al/(A0^Al*XX)  ) 

DRl  *  RlX*»(2./XXvl./Xl^2./XM^Bl/(B0^Rl*XX)  ) 

DO  =  PX«(2./XX42./X1>2./XM) 

DO  *  QX« (2./XX>2./Xlf (ni^2.*D2*XX)/(D0^Dl*XX*D2*XX*XX) ) 
DY20X  =  2.*RO«»ORO  +  2.*R1*DR1*CP*DP^DO 
DYDXd)  =  .5*DY2DX/(Y  d)*ELD*£LD) 

DYDX(NPTS)  =  OYOX(NPTS-l) 

DS  =  DX*SQRT(1.^0YDXd)*DYDXd)  ) 

IF  d.NE.2)  GO  TO  19 

5(2)  =  SORT  (  (X  (2) -X  (1 )  )  *(X  (2) -X(l)  )  ♦  (Y  (2) -Y  (1 )  )  *  (Y  (2)-Y(  1 )  )  ) 
GO  TO  27 
19  continue 

Sd)  =  S(I-l)  ♦  OS 
27  continue 
ST  =  S(I) 

TPd)  =  OYOXd)/(Y(I)*ELD) 
at  =  AT42.*PI*Yd)*0S 

2  continue 

WRITE(6.6)(  I*  Xd).  Yd)»  DYOX  ( I )  tS  ( I )  « TP  (I )  d  =  1  ♦NPTS ) 

G  F0RMAT(2X»  I3»  5F15.5) 

WRlTE(bO)  (Xd)  tYd)  f  I  =  1»NPTS) 

ENDFILE  60 

WRITE  (61)  (OYOXd)  ♦I  =  1,NPTS) 

VL3  =  .25<»PI*CP/(EL0*ELD) 

WRITE(6,15)  VL3 

15  F0RMAT(/1X.*  V0LUME/L3  =  **F1B.7/) 

WRlTE(8a)  VL3»AT»  ELO»  NPTS»  NCPTS 
WRITE  (6f25)  AT 

25  F0RMAT(/1X**  WETTED  AREA/L2  =  **  FlO.6) 

WRITE  (6»4)  ST 

4  format  (//»  1X»  *  TOTAL  ARC  LENGTH/L  =  *  ♦  F10.6) 

G0T07 

8  continue 

END 
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APPENDIX  C 


DESCRIPTION,  GLOSSARY,  LISTING,  AND  SAMPLE 
RUN  OF  DPIN3 


DESCRIPTION  OF  DPIN3 

20 

This  program  uses  the  polynomials  of  Granville  to  compute  the 
offsets  and  slopes  of  a  streamlined  body  of  revolution  formed  by  a  fore¬ 
body  with  rounded  nose,  a  parallel  mlddlebody  of  maximum  diameter,  and 
an  afterbody  with  pointed  tall.  The  axis  of  symmetry  is  the  x-axis,  the 
nose  is  at  the  origin,  and  the  tail  is  at  x  =  1.0.  On  the  forebody,  y 
is  equal  to  the  square  root  of  a  quintlc  in  x;  on  the  middlebody,  it  is 
(by  definition)  constant;  and  on  the  afterbody,  it  is  equal  to  the 
square  root  of  a  septic.  Coefficients  in  the  quintic  and  septic  are 

specified  by  geometric  parameters  which  are  input  variables,  together 

2  2 

with  the  requirements  that  the  slope  dy/dx  and  the  curvature  d  y/dx  be 
continuous  at  the  forebody-middlebody  and  middlebody-af terbody  junc¬ 
tures,  and  hence  be  zero  there. 

As  calculated  by  this  program,  the  bodies  are  represented  by  series 
of  frustums  of  cones  with  generating  axis  along  the  axis  of  revolution, 
since  the  body  shapes  are  formed  by  straight-line  segments  between 
points  defined  by  pairs  of  x-  and  y-  coordinates.  One  feature  of  the 
program  is  that  the  segment  lengths  or  step  sizes  over  the  foremost 
portion  of  the  body  are  one-tenth  as  large  as  those  over  the  rest  of  it. 
(A  different  ratio  of  step  sizes  could  be  attained  by  a  minor  modifica¬ 
tion  to  the  appropriate  DPIN  and  to  DPOUT.)  The  change  from  the  smaller 
to  the  larger  step  size  is  made  over  a  total  of  three  steps  in  such  a 
way  that  no  step  is  more  than  twice  as  large  as  the  one  preceding  it. 
Experience  in  operation  of  the  axisymmetric  potential-flow  program  has 
shown  the  desirability  of  this  relatively  gradual  change. 

Input  variables  to  the  program  consist  of  the  integers  NPTS,  NCPTS, 
and  INFL  and  the  floating-point  numbers,  MF,  MA,  KfL,  KAL,  RL,  SL,  and 


20 

Granville,  P,  S.,  "Geometrical  Characteristics  of  Noses  and  Tails 
for  Parallel  Middle  Bodies,"  NSRDC  Report  3763  (Dec  1972);  also  Inter¬ 
national  Shipbuilding  Progress,  Vol.  21,  No.  233,  pp.  3-19  (Jan  197A). 
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ELD.  NPTS  and  NCPTS  are  the  total  number  of  points  on  the  body  and  the 
number  of  closely  spaced  points.  INFL  is  the  maximum  number  of  inflec¬ 
tion  points  permitted  on  the  afterbody  and  must  be  either  0  or  1.  MF  is 
Che  ratio  of  forebody  length  to  total  length,  KFL  is  the  rate  of  change 
(with  x)  of  curvature  at  the  forebody-middlebody  juncture,  KAL  is  the 
rate  of  change  of  curvature  at  the  raiddlebody-af terbody  juncture,  RL  is 
the  radius  of  curvature  at  the  nose,  SL  is  the  slope  at  the  tail,  and 
ELD  is  the  ratio  of  total  length  to  maximum  diameter.  KFL,  KAL,  RL,  and 

SL  are  expressed  in  coordinates  scaled  so  that  the  total  length  of  the 

20 

body  is  1.  In  order  to  use  the  polynomials,  they  are  transformed  to 

stretched  coordinates  in  which  the  forebody  or  afterbody  length  and  the 

maximum  r.idius  each  equals  1.  Thus  the  coefficients  in  the  quintic  for 

the  forebody  are  R,  the  stretched  nose  radius,  and  KF,  the  stretched 

rate  of  change  of  curvature  at  the  forebody-middlebody  juncture.  The 

coefficients  in  the  septic  in  the  afterbody  are  KA,  the  stretched  rate 

of  change  of  curvature  at  the  middlebody-af terbody  juncture,  and  S2,  the 

square  of  SP,  the  stretched  slope  at  the  tail.  KA  and  SP  have  their 

20 

signs  reversed  so  as  to  be  ordinarily  positive.  Granville  avoids 

undesirable  bulges,  zeros,  and  inflection  points  on  the  forebody  by 

requiring  R  and  KF  to  lie  in  an  admissible  region  when  graphed  (see 

Figures  2-4  of  Reference  20,  in  which  they  are  denoted  by  r  and  k^^). 

In  the  same  way,  undesirable  afterbodies  are  excluded  by  requiring  that 

S2  and  KA  lie  in  one  of  two  admissible  regions  when  graphed  (see  Figures 

2 

5  and  6  of  Reference  20,  in  which  they  are  denoted  by  s  and  k^),  accord¬ 
ing  to  whether  zero  or  one  inflection  points  are  to  be  allowed  on  the 
afterbody.  Instead  of  introducing  the  complications  associated  with 
numerical  representation  of  the  curves  which  delineate  the  admissible 
regions  on  these  graphs,  DPIN3  uses  a  different  method  of  checking  the 
occurrence  of  bulges  and  points  of  Inflection:  as  the  slope  at  each 
point  is  computed,  it  is  compared  with  zero  and  the  slope  at  the  pre¬ 
vious  point.  If  at  points  on  the  forebody  the  slope  is  negative  or 
greater  than  the  previous  slope,  the  program  prints  out  a  message  that 
the  forebody  parameters  are  inadmissible  and  stops.  At  points  on  the 
afterbody,  the  program  prints  out  the  corresponding  message  and  stops 
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if  the  slope  is  positive  or,  if  no  inflection  points  are  permitted  (INFL 
■  0),  the  slope  is  greater  than  the  previous  slope.  If  one  inflection 
point  is  permitted  (INFL  ■  1),  the  slope  is  allowed  to  exceed  the  prev¬ 
ious  slope  over  a  range  of  points,  past  which  it  may  become  less  than 
the  previous  slope  over  one  range  but  must  not  exceed  the  previous  slope 
over  a  second  range. 

The  program  user  may  wish  to  modify  DP1N3  so  that  KF,  KA,  R,  and  S2 

(or  its  square  root)  are  used  as  input  rather  than  the  unstretched  KKl., 

KAL,  RL,  and  SL.  Furthermore,  he  may  wish  to  replace  either  KK  or  R 

with  CDF,  the  forebody  prismatic  coefficient,  or  to  replace  either  KA  or 

S2  with  CPA,  the  afterbody  prismatic  coefficient  since  Equations  (103) 

20 

and  (104)  in  Granville  present  linear  relations  between  these  re¬ 
spective  triads  of  geometric  parameters. 

After  DPIN3  has  read  its  ten  input  variables,  it  writes  them  out, 
then  computes  the  forebody,  afterbody,  and  overall  prismatic  coefficients 
and  writes  them  out.  Then  it  finds  the  appropriate  step  size  DX  and 
proceeds  to  calculate  the  offsets  Y(I)  at  each  point  along  the  forebody 
and  afterbody  by  taking  the  square  root  of  the  appropriate  polynomial. 

The  slope  DYDX(I)  and  increment  in  arc  length  DS  are  calculated  by  using 
the  derivative  of  the  polynomial.  The  wetted  area  and  a  geometric  func¬ 
tion  TP(1),  which  is  used  for  predicting  boundary-layer  transition,  are 
also  computed  at  each  point.  Then  I,  X,  Y,  DYDX,  S,  and  TP  are  written 
out  for  each  point.  Finally  the  volume,  and  forebody,  middlebody, 
afterbody,  and  overall  wetted  areas  and  arc  lengths  are  written  out. 

X(I),  Y(I),  and  the  length-to-diameter  ratio,  and  numbers  of  total 
points  and  closely  spaced  points  are  also  written  on  tape  for  use  by 
subsequent  programs. 

DPIN3  -  GLOSSARY 

The  following  glossary  of  variables  used  in  DPIN3  is  arranged 
alphabetically  by  FORTRAN  variable  name. 


FORTRAN 

Variable 

Name 

Variable 

Def inition 

AA 

Wetted  area  of  afterbody  divided  by 
square  of  total  length 

AF 

a<«/l2 

Wetted  area  of  forebody  divided  by 
square  of  total  length 

AM 

Wetted  area  of  middlebody  divided  by 
square  of  total  length 

AT 

Total  wetted  area  divided  by  square  of 
total  length 

CP 

c 

p 

Prismatic  coefficient  of  complete  body 

CPA 

c 

p 

Prismatic  coefficient  of  afterbody 

CPF 

c 

p 

Prismatic  coefficient  of  forebody 

DS 

ds/L 

Increment  in  arc  length  s  divided  by 
total  length 

DX 

dx/L 

Increment  in  axial  length  x  divided  by 
total  length 

DYDX 

dx 

Slope;  dimensioned  to  be  a  function  of 
I 

ELD 

L/MAX.  DIAMETER 

Ratio  of  total  length  to  maximum 
diameter 

I 

i 

Integer  which  increases  from  1  at  the 
nose  to  NCPTS  at  the  point  where  step 
size  starts  increasing,  hence  to  NPTS 
at  the  tail 

INFL 

Control  variable  which  specifies  the 
maximum  number  of  Inflection  points 
allowed  on  the  afterbody;  always  input 
as  0  or  1 
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FORTRAN 

Variable 

Name 

Variable 

Def Inlt Ion 

KA 

kj^(a) 

Rate  of  change  of  curvature  at 
mlddlebody-af terbody  juncture,  normal¬ 
ized  for  unit  afterbody  length  and 
unit  maximum  radius,  with  sign 
reversed  so  as  to  be  ordinarily 
positive 

KAL 

dx' 

(mlddlebody- 
af terbody 
juncture) 

Rate  of  change  of  curvature  at 
middlebody-af terbody  juncture  times 
square  of  total  length 

KF 

k 

‘"l 

Rate  of  change  of  curvature  at 
forebody-middlebody  juncture,  normal¬ 
ized  for  unit  forebody  length  and  unit 
maximum  radius 

KFL 

2  di 

dx^ 

(forebody- 

middlebody 

juncture) 

Rate  of  change  of  curvature  at 
forebody-middlebody  juncture  times 
square  of  total  length 

MA 

Ratio  of  length  of  forebody  plus 
middlebody  to  total  length 

MAI 

1 

1  -  m 

Ratio  of  afterbody  length  to  total 
length 

MF 

Ratio  of  forebody  length  to  total 
length 

NCPTS 

Number  of  the  point  downstream  of 
which  the  point  spacing  begins  to 
Increase 

NCI 

-b 

1 

NC2 

2 
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FORTRAN 

Variable 

Name 

NC3 

NC4 

NPTS 

PI 

R 

RL 

S 

SA 

SF 

SL 

SM 

SP 

ST 


Variable 


Def Initlon 


+  3 


4-  A 

n  +4 


1 

L 


(nose) 


s/L 


s(«/L 


dx 


(tail) 


8<"'/L 


,(P) 


s 


(t) 


Total  number  of  points 

Ratio  of  the  circumference  of  a  circle 
to  its  diameter 

Radius  of  curvature  at  nose,  normal¬ 
ized  for  unit  forebody  length  and  unit 
maximum  radius 


Radius  of  curvature  of  nose  divided  by 
total  length 


Arc  length  at  arbitrary  x  divided  by 
total  length;  dimensioned  to  be  a  func¬ 
tion  of  I 


Arc  length  of  afterbody  divided  by 
total  length 


Arc  length  of  forebody  divided  by 
total  length 


Slope  at  tail 


Arc  length  of  middlebody  divided  by 
total  length 


Slope  at  tail,  normalized  for  unit 
afterbody  length  and  unit  maximum 
radius,  with  sign  reversed  so  as 
to  be  ordinarily  positive 

Total  arc  length  divided  by  total 
length 
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FORTRAN 

Variable 

Name 

S2 

TP 

VL3 

X 

XA 

XAl 

XF 

XFl 

XX 

Y 


Variable 


Definition 


,(P)2 


SP  squared 


1  dy 
ELD  y(x)  dx 


Function  used  in  DPOUT  to  predict 
transition;  dimensioned  to  be  a  func¬ 
tion  of  I 


Volume  of  body  divided  by  cube  of 
total  length 


x/L 


Axial  coordinate  divided  by  total 
length;  dimensioned  to  be  a  function 
of  I 


Reversed  x,  normalized  to  increase 
from  XA  =  0  at  tail  to  XA  -  1  at 
middlebody-af terbody  juncture 


x/L 


.(a) 


Reversed  x,  normalized  to  increase 
from  XAl  =  -1  at  tail  to  XAl  =  0  at 
middlebody-af terbody  juncture 


x/L 

_(f) 


X  normalized  to  increase  from  XF  =  0 
at  nose  to  XF  =  1  at  forebody- 
middlebody  juncture 


x/L 


m 


(f) 


m 


(f) 


y/L 


X  normalized  to  increase  from  XFl  =  -1 
at  nose  to  XFl  =  0  at  forebody- 
middlebody  juncture 

Axial  variable,  indentical  to  X  except 
that  it  is  not  dimensioned  as  a 
function  of  I 


Body  radius  divided  by  total  length; 
dimensioned  to  be  a  function  of  1 
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LISTING  AND  SAMPLE  RUN  OF  DPIN3 
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APPENDIX  D 


DESCRIPTION  OF  DA50  AND  DESCRIPTION,  GLOSSARY,  LISTING, 
AND  SAMPLE  RUN  OF  DPOUT 


DESCRIPTION  OF  DA50 

DA50  is  the  DTNSRDC  designation  of  the  Douglas-Neumann  potential- 

12 

flow  program  for  axisymmetr ic  bodies.  This  program  represents  bodies 
by  ensembles  of  flat  quadrilaterals,  on  each  of  which  there  is  an 
(initially)  unknown  constant  source  strength.  A  Fredholm  integral  equa¬ 
tion  of  the  second  kind  for  the  velocity  potential  results  from  requiring 

the  flow  to  be  tangent  to  the  body  surface.  (This  is  Equation  (2.3)  in 
12 

Hess  and  Smith.  )  Since  there  are  only  a  finite  number  of  panels  and 
sources  which  Induce  the  velocity  potential,  the  integral  equation  is 
replaced  by  a  system  of  linear  algebraic  equations,  of  number  equal  to 
the  number  of  quadrilaterals,  each  one  of  which  involves  all  of  the 
source  strengths.  Thus  solution  of  the  algebraic  equations  involves 
inversions  of  square  matrices  with  the  lengtn  of  each  side  equal  to  tlie 
number  of  quadrilaterals.  After  the  source  strengths  are  determined, 
the  velocity  potential  and  velocity  components  and  pressure  anywhere  in 
the  flow  field  can  easily  be  found  in  terms  of  the  source  strengths. 

In  the  axisymmetric  case,  the  quadrilaterals  become  frustums  of 
cones,  with  axes  along  the  body's  axis  of  symmetry.  DA50  is  so  long 
that  it  is  divided  into  six  overlay  links.  LINK  1  sets  up  the  panels 
and  algebraic  equations  tor  the  source  strengths.  LINK  2,  LINK  3, 

LINK  4  and  LINK  5  provide  options  for  solving  the  algebraic  equations. 

LINK  2,  which  uses  the  Seidel  iteration,  is  the  only  one  which  has  been 
used  in  this  project.  (LINK  3  solves  the  algebraic  equations  by  a 
direct  method;  LINK  4  prepares  a  tape  for  LINK  5  which  solves  the  matrices 
by  successive  orthogonalization. )  LINK  b  computes  the  desired  velocity 
potential,  velocity  components,  and  pressure  coefficients  due  to  the 
source  strengths  and  prints  them.  Thus  when  DA50  is  entered  again  to 
compute  the  revised  velocity  distribution  due  to  the  effect  of  displace¬ 
ment  thickness,  only  LINK  6  is  used;  no  repetition  of  the  time-consuming 
inversion  is  needed. 
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DESCRIPTION  OF  DPOUT 


This  program  uses  Che  results  of  DPIN  and  DA50  together  with  viscous- 
flow  theory  presented  elsewhere  in  this  report  and  in  References  1-4  to 
compute  the  boundary  layer  on  a  body  of  revolutioi;.  The  boundary  layer 
is  laminar  at  the  nose,  undergoes  transition  to  turbulence  at  some  point 
downstream  of  the  nose  if  the  overall  Reynolds  number  R^^  is  sufficiently 
large,  continues  as  a  turbulent  boundary  layer  to  the  tall,  and  leaves 
the  body  as  a  turbulent  wake.  However,  if  R^^  is  sufficiently  small  that 
laminar  separation  occurs  on  the  forebody,  transition  can  be  forced  to 
occur  at  the  lamlnary-separ:  tion  point.  DPOUT  includes  four  different 
methods  of  controlling  transition: 

1.  The  method  of  Granville^  is  generalized  to  lower  values  of 

"o  d 

- r^,  where  D  is  the  maximum  diameter,  and  to  flows  with  high  back- 

y  dx  o  ® 

ground  turbulence.  This  includes  three  subcases. 

1  3 

2.  The  method  of  Granville  as  used  by  Smith. 

3.  Transition  can  be  forced  to  occur  at  a  designated  point,  either 
because  it  has  been  observed  to  occur  there  in  an  experiment  or  because 
a  trip  has  been  placed  there.  If  a  trip  is  used,  DPOUT  includes  an 
option  increase  in  R^^  (and  ultimately  in  total  drag)  to  account  for  the 
parasitic  drag  of  the  trip.  Designation  of  a  transition  point  does  not 
prevent  transition  from  occurring  upstream  of  it  if  this  is  predicted  by 
one  of  the  other  three  options  which  has  been  chosen. 

4.  If  laminar  separation  occurs  upstream  of  the  transition  point 
determined  by  one  of  the  other  three  methods,  the  user  has  the  option  of 
defining  the  separation  point  as  the  transition  point.  Thus  when  this 
option  is  selected,  it  is  assumed  that  a  separated  laminary  boundary 
layer  immediately  reattaches  as  a  turbulent  boundary  layer. 

The  program  also  includes  the  option  of  setting  the  velocity  and 
wall  radius  at  constant  values  and  thus  calculating  the  flat-plate  drag 
as  a  check.  In  this  case  the  body  is  effectively  replaced  by  a  circular 
cylinder  with  a  very  thin  wall  and  a  radius  equal  to  the  maximum  body 
radius,  and  the  drag  on  the  inner  surface  of  the  wall  is  not  included. 


Input  variables  to  DPOUT  include,  aside  from  the  data  on  tape 

written  by  DPIN  and  DA50,  the  following: 

ICONTRL  This  integer  a  control  variable  which  should  be  zero  the 
first  time  the  boundary  layer  is  calculated  for  a  given  body 
and  overall  length  Reynolds  number  R^,  if  a  recomputation  is 

to  be  made  which  includes  the  augmented  hydrodynamic  source 
strength  to  represent  the  effect  of  the  displacement  thick¬ 
ness.  If  no  recomputation  is  to  be  made,  set  ICONTRL  =  1. 

If  the  run  is  to  use  the  recomputed  data  (i.e.  it  is  not 
an  initial  run  to  obtain  DELSIG),  also  set  ICONTRL  =  1. 

153  This  integer  is  a  control  variable  which  should  be  1  if  and 

only  if  it  is  desired  to  use  the  transition  criterion  of 
Reference  1  Instead  of  that  of  Reference  2. 

ITRIP  This  integer  is  a  control  variable  which  should  be  1  if  and 
only  it  J.t  is  desired  that  transition  take  place  at  XTRIP  if 
it  has  not  occurred  upstream  of  it. 

XTRIP  This  floating-point  number  is  the  location  where  transition  is 

forced,  expressed  as  axial  distance  from  the  nose  divided  by  total 
length.  It  must  be  input  even  if  ITRIP  is  unequal  to  1. 

IDRTH  This  integer  is  a  control  variable  which  should  be  1  if  and 
only  if  it  is  desired  that  at  the  point  where  transition  is 
forced,  an  increment  be  added  to  the  momentum-thickness 
Reynolds  number  to  represent  the  added  drag  due  to  a  sand 
strip.  It  must  be  input  even  if  ITRIP  is  unequal  to  1, 

IFP  This  integer  is  a  control  variable  which  should  be  1  if  and 

only  if  it  is  desired  to  replace  the  body  shape  by  a  circular 
cylinder  with  its  surface  parallel  to  the  undisturbed  flow 
and  to  replace  the  computed  velocity  by  the  undisturbeH-f low 
speed.  This  provides  flat-plate  results  which  give  a  useful 
check  on  the  computations. 

ITLS  This  integer  is  a  control  variable  which  should  be  1  if  and 

only  if  it  is  desired  to  force  transition  to  take  place  at  the 
laminar  separation  point  in  cases  where  the  Reynolds  number  is 
so  low  that  transition  would  not  otherwise  occur  until  down¬ 
stream  of  the  laminar  separation  point.  Thus  it  is  assumed 
that  a  separated  laminar  boundary  layer  immediately  reattaches 
as  a  turbulent  boundary  layer. 

RL  This  floating-point  number  is  the  Reynolds  number  based  on  bodv 

length  and  velocity  of  the  undisturbed  flow. 


J  This  integer  is  a  control  variable  which  is  1  if  and  only  if, 

when  using  the  transition  criterion  of  Reference  2,  it  is 
desired  to  use  the  low-background-turbulence  correlation  curve. 

It  should  be  input  even  if  153  is  1. 

L  This  Integer  is  a  control  variable  which  is  1  if  and  only  if, 

when  the  transition  criterion  of  Reference  2  is  used  with 
J  =  1  and  TP  is  sufficiently  large,  it  is  desired  to  use  as 
the  transition  criterion  a  constant  difference  of  639  between 
RTHT  and  RTHN.  Although  it  is  meaningless  with  J  not  equal  to  1 
or  if  153  =  1,  it  should  be  input  in  these  cases. 

The  program  is  able  to  treat  the  same  body  with  a  number  of  differ¬ 
ent  values  of  RL,  J,  and  L  because  after  execution  for  one  set  of  values, 
it  returns  and  reads  the  next. 

As  the  program  reads  NPTS,  NCPTS,  VL3,  AT,  153,  ITRIP,  XTRIP,  IDRTH, 
IFP,  and  ITLS,  it  writes  them  out.  It  also  reads  the  body-shape  descrip¬ 
tion  calculated  by  DPIN  and  the  velocity  distribution  calculated  by  DA50. 
DA50  gives  the  tangential  velocity  at  the  center  of  the  integration  steps, 
and  DPOUT  uses  velocity  at  the  ends  of  the  steps,  so  the  latter  is  found  by 
linear  interpolation.  Linear  interpolation  of  velocity  at  95  percent  of 
the  axial  length  is  also  used  over  the  last  5  percent  of  the  body.  The 
x-derivative  of  velocity  is  also  needed  by  the  program,  and  this  is  calcu¬ 
lated  by  a  numerical-differentiation  subroutine,  DGT3,  which  is  described 
in  the  IBM  reference  manual.  If  IFP  is  1,  Y(I)  is  replaced  by  the  maxi¬ 
mum  radius  and  U(I)  by  the  speed  of  the  undisturbed  flow.  Next  I,  X,  Y, 

U,  and  DUDX  are  written  out  for  points  on  the  body.  RL,  J,  and  L  are 
read  for  the  case  and  written  out.  Various  quantities  used  in  computing 
the  boundary  layer  are  initialized  and  the  laminar-boundary-layer  calcula¬ 
tions  are  begun  by  entering  a  DO  loop  in  which  the  neutral-stability  point 
is  located  if  RL  is  large  enough  for  there  to  be  such  a  point.  The 
methods  by  which  the  laminar  boundary  layer  and  the  neutral-stability 
point  are  calculated  are  described  elsewhere.  If  no  neuti a] -stability 
point  is  found,  a  message  to  the  effect  that  RL  is  too  small  is  printed 
out  and  the  program  goes  to  the  next  RL,  J,  and  L.  If  the  neutral- 
stability  point  is  found,  the  program  writes  out  I,  X,  Y,  U,  RTH,  etc.  at 
that  point,  denoting  them  with  N  as  the  last  letter,  and  continues  the 


90 


laminar-boundary-iayer  calculations  by  entering  a  new  DO  loop  In  which 
the  transition  point  is  located  by  the  method  of  Reference  2  or,  if  153 
is  1,  by  the  method  of  Reference  1  as  long  as  RL  la  large  enough  for 
there  to  be  such  a  point.  If  RL  Is  not  large  enough,  the  program  prints 
out  a  message  to  that  effect  and  goes  to  the  next  RL,  J,  and  K.  Details 
of  the  two  transition-prediction  methods  have  been  summarized  previously. 
Both  upstream  and  downstream  of  the  neutral-stability  point  the  shear 
stress  is  calculated  according  to  the  Blasius  flat-plate  formula,  and  if 
ITRIP  is  1,  transition  takes  place  immediately  If  X  passes  XTRIP.  Also, 
both  upstream  and  downstream  of  the  neutral  stability  point  a  check  for 
laminar  separation  is  made.  If  laminar  separation  is  encountered,  the 
program  prints  out  a  message  to  chat  effect  and  goe^  to  the  next  RL,  J, 
and  K  unless  ITLS  is  1,  in  which  case  transition  takes  place  at  the 
laminar-separation  point.  I,  X,  Y,  U,  PGP,  and  RTH  are  written  out  at 
each  point  on  the  laminar  boundary  layer. 

At  the  transition  point,  I,  X,  Y,  U,  RTH,  etc,  are  again  printed 

out,  denoted  by  T  as  the  last  letter,  and  the  turbulent-boundary-layer 

calculations  are  begun  by  entering  a  new  DO  loop.  H  at  transition 

(denoted  by  HT)  is  given  by  Equation  (15).  The  method  by  which  the 

turbulent  boundary  layer  is  calculated  appears  elsewhere  in  this  report 

4 

and  has  been  given  in  detail  by  Granville.  The  calculations  begin  at 
the  first  point  downstream  of  the  transition  point  through  which  it  is 
assumed  that  TH,  and  hence  RTH,  are  continuous.  Then  HT  is  found  by 
using  Equation  (15)  and  SG  is  found  by  using  Equation  (16).  The  method 
by  which  the  transcendental  Equation  (16)  is  solved  for  SG  with  H  and 
RTH  known  consists  of  using  H  and  RTH  to  evaluate  its  right-hand  side 
(denoted  by  FR) .  As  a  first  try,  SG  =  11.0  is  used.  If  the  left-hand 
side  (denoted  by  FS)  is  less  than  FR,  SG  is  increased  by  100  and  FS  is 
reevaluated.  If  FS  is  now  no  less  than  FR,  SG  is  decreased  in  steps  of 
10.0  until  it  is  less  than  FR.  In  the  same  way,  if  the  FS  based  on 
SG  =  11.0  had  been  no  less  than  FR,  SG  would  be  decreased  by  steps  of 
10.0  until  FS  were  less  than  FR.  The  process  is  continued  in  this  way; 

SG  increases  in  steps  of  100,  decreases  in  steps  of  10,  increases  in 
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steps  of  1,  etc.,  with  comparisons  of  FS  and  FR  made  at  each  stage, 
until  SG  has  been  found  to  an  accuracy  of  0.001.  Next  Equations  (17)- 

(34)  are  solved  so  that  the  right-hand  sides  of  Equations  (35)  and  (36) 
can  be  evaluated.  Equations  (35)-(38)  are  solved  by  using  a  two-stage 
predictor-corrector  technique,  governed  by  a  control  variable  K  which  is 
1  in  the  first  ("predictor")  stage  and  2  is  the  second  ("corrector") 
stage.  With  K  =  1  the  coefficients  and  right-hand  sides  of  Equations 

(35)  through  (38)  are  evaluated  by  using  quantities  at  the  left-hand 
ends  of  the  integration  steps.  The  coefficients  and  Increments  are 
denoted  by  D0ML2P,  DPSIL2P,  AlP,  BIP,  CIP,  A2P,  B2P,  C2P,  C3P,  C4P, 

DTHP,  and  DHP.  After  they  have  been  evaluated  and  0ML2,  PSIL2,  TH,  and 
H  incremented  by  D0ML2P,  DPSIL2P,  DTHP,  and  DHP,  the  program  changes  K 
to  2  and  returns  to  the  point  where  SG  is  found  from  Equation  (16)  and 
repeats  the  calculations.  During  this  stage,  the  right-hand  sides  of 
Equations  (35)-(38)  are  evaluated  by  using  values  of  0ML2,  PSIL2,  TH, 

and  H  which  were  estimated  at  the  "predictor"  stage.  The  revised  coeffi¬ 
cients  and  increments  are  denoted  by  D0ML2C,  DPSIL2C,  AlC,  BIC,  CIC, 

A2C,  B2C,  C2C,  C3C,  C4C,  DTHC,  and  DHC.  The  final  increments  in  0ML2, 
PSIL2,  TH  and  H  over  the  values  at  the  beginning  of  the  integration  step 
are  0.5  *  (D0ML2P  +  D0ML2C),  0.5  *  (DPSIL2P  +  DPSIL2C),  0.5  *  (DTHP  + 

DTHC),  and  0.5  *  (DHP  +  DHC).  The  shear-stress  coefficient  CF  is  calcu- 

2 

lated  by  using  Equation  (4),  with  C^  now  equal  to  2/SG  .  After  all  the 
calculations  have  been  made  at  this  point  on  the  turbulent  boundary 
layer,  the  values  of  I,  0ML2,  PSIL2,  ELSL2,  TH,  H,  SG,  log^  (RTH) ,  and 
CF  are  written  out.  Next  the  program  returns  to  the  beginning  of  the  DO 
loop,  sets  K  back  to  1,  and  repeats  the  calculations  on  the  second 
integration  step  downstream  of  transition.  The  only  difference  is  that 
0ML2  and  PS1L2  at  the  beginning  of  the  new  step  are  now  known  and  Equa¬ 
tions  (32)  and  (33)  are  not  used.  The  process  is  continued  until  the  DO 
loop  is  completed  and  quantities  have  been  computed  at  each  point  along 
the  turbulent  boundary  layer.  However,  a  different  computation  method 
is  used  over  the  last  five  points;  0ML2,  PSIL2,  TH,  and  H  are  computed 
by  linear  extrapolation  of  the  values  at  the  two  previous  points. 
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After  completion  of  the  DO  loop  for  the  turbulent  boundary  layer, 
DELSIG  is  computed  at  each  point  along  the  body  if  ICONTRL  is  not  1. 
DELSIG  is  computed  by  using  Equation  (6)  except  on  the  last  five  points, 
where  linear  interpolation  from  the  two  previous  points  is  used. 

The  values  of  quantities  at  the  tail,  namely,  1,  HE,  THE,  0ML2E, 

FCT  2F 

PS1L2E,  ELSL2E,  SHE  (defined  by  SHE  *  ;rr) ,  and  UE  are  written  out. 

UML2h 

0ML2D  is  computed  by  using  Equations  (48),  (49),  and  (50).  The  drag 

2/3 

coefficients  based  on  (volume)  and  wetted  area  as  reference  areas  are 
computed  by  using  Equation  (47)  and  written  out.  If  ICONTRL  Is  not  1, 
values  of  DELSIG  at  the  centers  at  the  integration  steps  are  computed  by 
using  linear  interpolation  from  the  end  points  of  the  steps.  These  are 
stored  on  tape  for  use  during  the  reexecution  of  last  part  of  DA50  to 
represent  the  effect  of  displacement  thickness. 

DPOUT  -  GLOSSARY 

The  following  glossary  of  variables  used  in  DPOUT  is  arranged 
alphabetically  by  FORTRAN  variable  name. 


FORTRAN 

Variable 

Variable 

Definition 

Name 

ALOGPGF 

anti  logjQ 

(PGF) 

Inverse 

of 

logjQ  (PGF) 

ALRTH 

anti  log 

e 

(RTH) 

Inverse 

of 

log^  (RTH) 

AT 

A<^) 

Total  wetted  area 

AlC 

A 

^1 

Element 

in 

matrix  for  DHC 

and 

DTHC 

AlP 

A 

^1 

Element 

in 

matrix  for  DHP 

and 

DTHP 

A2C 

A 

Element 

in 

matrix  for  DHC 

and 

DTHC 

A2P 

A 

Element 

in 

matrix  for  DHP 

and 

DTHP 
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FORTRAN 

Variable  Variable 
Name 


B 

6 

BE 

g(e) 

BLAM 

X 

BO 

e 

o 

BIG 

BIP 

B^^P) 

B2C 

B 

®2 

B2P 

CD 

S 

CF 

CFN 


9A 


fortran' 


Variable 

Name 

Variable 

Def inlt ion 

CFT 

.  (t) 

at  transition 

CS 

C 

s 

Drag  coefficient  based  on  wetted  area 

CIC 

P  (c) 

^1 

(  c  ) 

Term  used  in  calculating 

CIP 

.  (P) 

^1 

Term  used  in  calculating 

C2C 

^  (c) 

^2 

Term  used  in  calculating 

C2P 

c  (P) 

2 

Term  used  in  calculating 

C3C 

c 

^3 

Element  in  matrix  for  DHC  and  DTHC 

C3P 

P  (P) 

^3 

Element  in  matrix  for  DHP  and  DTHP 

CAC 

c 

Element  in  matrix  for  DHC  and  DTHC 

CAP 

c  (P) 

Element  in  matrix  for  DHP  and  DTHP 

x/L 

D 

D 

r  [y  (x'/L)]^  [sec  a  (x'/L)]  d  x/L 

X  /L 
n 

(evaluated  for  x  <  x  <  x„) 
n  —  —  1 

DCOS 

d  cos  a 

Increment  on  cos  rx  (x/L)  over  an 
integration  step 

DDELSIG 

DELSIG  evaluated  at  center  of  integration 
step  by  linear  interpolation;  dimensioned 
to  be  a  function  of  I 

DELDX 

dELSL2 

dx 

Derivative  of  ELSL2  with  respect  to  x; 
dimensioned  to  be  a  function  of  1 
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FORTRAN 

Variable 

Name 

DELSIC; 

DHC 

DHDH 

DHDR 

DHP 

DHPH 

DHPR 

DHR 

DHTH 

DOML2C 

DOML2P 

DPSIL2C 


Variable 


Def  Inlt  l.on 


6  SIC 


6  H 


(c) 


9  HU 
9  H 


Increment  In  source  strength  SIG  (of 
DA50)  which  represents  the  effect  of 
displacement  thickness;  dimensioned  to 
be  a  function  of  I 


Increment  in  H  over  an  Integration  step, 
evaluated  at  corrector  stage  of  predictor- 
corrector  method 


Derivative  of  HD  with  respect  to  H 


9  HD 

9  (log^  RTH) 


Derivative  of  HD  with  respect  to  log^  RTH 


6  H 


(p) 


9  HP 
9  H 


Increment  in  H  over  an  integration  step, 
evaluated  at  predictor  stage  of  predictor- 
corrector  method 


Derivative  of  HP  with  respect  to  H 


9  HP 

9  (log^  RTH) 


Derivative  of  HP  with  respect  to  log^  RTH 


9  HT 

9  (log^  RTH) 

9  HT 
9  H 


Derivative  of  HTL  with  respect  to  log^  RTH 
Derivative  of  HTL  with  respect  to  H 


6  0ML2^‘^^ 

6  0ML2^P^ 

6  PSIL2^‘^^ 


Increment  in  0ML2  over  an  integration 
step,  evaluated  at  corrector  stage  of 
predictor-corrector  method 


Increment  in  0ML2  over  an  integration 
step,  evaluated  at  predictor  stage  of 
predictor-corrector  method 


Increment  in  PSIL2  over  an  integration 
step,  evaluated  at  corrector  stage  of 
predictor-corrector  method 
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FORTRAN 

Variable 

Name 


Variable 


Definition 


Dl’SILzP 

DRTH 

DSH 

DSR 

DTHC 

DTHP 

DU 

DUDX 

DX 

DY 

UYDX 

DO 

E 

EC 

EE 


(S  PSIL2^P^ 


6  RTH 

L-  9  SG 

SG  3  H 

1 _ 3  SG 

SG  3  (log^  RTH) 

6  TH^^^ 

6  TH^P^ 


6  U 

M 

dX 

dX 

6  Y 

dY 

dX 

D 

o 

E 

E 

e 


Increment  In  PSIL2  over  an  Integration 
step,  evaluated  at  predictor  stage  of 
predictor-corrector  method 

Increment  In  RTH  due  to  the  presence  of  a 
sand  strip  used  as  a  transition  trip 

1/SG  times  derivative  of  SG  with  respect 
to  H 


1/SG  times  derivative  of  SG  with  respect 

to  log  RTH 
e 

Increment  in  TH  over  an  integration  step, 
evaluated  at  corrector  stage  of  predictor 
corrector  method 


Increment  in  TH  over  an  integration  step, 
evaluated  at  predictor  stage  of  predictor 
corrector  method 


Increment  in  U  over  an  integration  step 


Derivative  of  U  with  respect  to  X;  dimen¬ 
sioned  to  be  a  function  of  I 

Length  of  an  integration  step 

Increment  in  Y  over  an  integration  step 


Derivative  of  Y  with  respect  to  X;  dimen¬ 
sioned  to  be  a  function  of  I 

Value  of  D  at  previous  integration  step 

2 

Entrainment  factor,  equal  to  EC/SG" 

Reduced  entrainment  factor 
Base  of  natural  logarithms 
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FORTRAN 

Variable 

Name 


Variable 


Def Initlon 


ELD 


ELSL2 


ELSL2E 

F 

FD 

FDO 

FR 


FRTHN 


FS 


FO 

G 

H 

HD 

HE 

HP 

HT 


L/MAX.  diameter  Ratio  of  body  length  to  maximum  diameter; 

here  the  body  length  Is  one  so  it  is  the 
inverse  of  the  diameter 


A*/L^ 


ELSL2^®^ 

F 

FD 


FD 


F 


(r) 


F(RTH^"^) 


,(s) 


F 

o 

G 


H 

"a 

„<e) 

„(0 


Displacement  area  divided  by  square  of 
body  length 

ELSL2  evaluated  at  tail 
Integrand  of  B 
Integrand  of  D 

Value  of  FD  at  previous  integration  step 

Function  of  log  RTH  and  H  used  to  com- 
e 

pute  SG 
RTH2RL^"^ 


(n))  ^ 


Function  of  SG  and  H  used  in  computing 
SG  by  taking  trial  values  of  SG  and 
comparing  FS  with  FR 

Value  of  F  at  previous  integration  step 


Rotta  shape  parameter,  equal  to  SG 


Shape  parameter,  equal  to  6*/0 
Quadratic  displacement-shape  parameter 


H  evaluated  at  tail 

Quadratic  momentum-shape  parameter 


H  evaluated  at  transition  point 


k 
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FORTRAN 

Variable 

Name 


Variable 


Def inition 


HTL 

H 

Entrainment-shape  parameter 

1 

i 

Integer  which  increases  from  1  at  the 
nose  to  NCPTS  at  the  point  where  tlie 
step  size  begins  increasing,  hence  to 

NPTS  at  the  tall 

IBID 

IBID 

Dummy  variable  for  I,  used  in  DO  loop 
when  DELSIG  is  computed  in  the  case  where 

RL  is  so  small  that  the  point  of  neutral 
stability  is  never  reached  on  the  body 

ICAN 

ICAN 

Dummy  variable  for  I,  used  in  DO  loop 
when  DELSIG  is  computed  in  the  case  where 
transition  to  turbulence  occurs  on  the 
body 

ICONTRL 

ICONTRL 

Control  variable  which  is  1  when  the 
execution  of  DPOUT  uses  the  velocity 
distribution  corresponding  to  the 
hydrodynamic  source  strength  which 
includes  the  effect  of  displacement 
thickness 

IDRTH 

^(dRTH) 

Control  variable  which  is  meaningful  only 
if  ITRIP  is  1;  IDRTH  is  1  if  and  only  if 
it  is  desired  to  increment  RTH  to  represent 
added  drag  due  to  the  effect  of  a  sand 
strip  which  stimulates  transition 

lER 

lER 

Error  parameter  used  in  subroutine  DGT 3 

IFF 

j(FP) 

Control  variable  which  is  1  if  and  only 
if  it  is  desired  to  obtain  flat-plate 
results  by  setting  U  =  1.0  and 

Y  =  (maximum  diamctor/2) 

IN 

j  (n) 

1  at  neutral  stability  point 

INI 

l(n)  ^  j 

FORTRAN 

Variable 

Name 

[QUIT 

IT 

ITLS 

ITRIP 

ITRY 

ITl 

153 

J 


Variable 


Def Inltion 


I  QUIT 


j(tls) 


, (trip) 


ITRY 


+  1 


Control  variable  which  ensures  termination 
and  procession  to  the  next  R^^  in  a  case 

wlu're  laminar  separation  occurs  for  an  Rj 


I  at  transition 


Control  variable  which  is  1  if  and  only 
If  it  is  desired  that  if  laminar  separa¬ 
tion  occurs  at  a  point,  then  that  point 
will  be  defined  as  the  transition  point 

Control  variable  which  is  1  if  and  only 
if  it  is  desired  to  force  transition  at  a 
point  (x  =  XTRIP)  if  it  has  not  occurred 
upstream  of  that  point 


Dununy  variable  for  1,  used  in  the  DO  loop 
where  DF>LSIG  is  computed  in  the  case  where 
RL  is  so  small  that  the  transition  point 
is  not  reached 


I  at  the  first  point  downstream  of  the 
transition  point 


(53) 


K 


Control  variable  which  should  be  set  equal 
to  1  if  and  only  if  it  is  desired  to  use 
the  transition  criterion  of  Reference  1 
rather  than  that  of  Reference  2 

Control  variable,  meaningful  only  when 
the  transition  criterion  of  Reference  2 
is  being  used;  J  =  1  gives  the  correlation 
curve  for  low-background  turbulence  and 
J  =  2  gives  the  correlation  curve  for  high- 
background  turbulence 

Control  variable  for  predictor-corrector 
method  of  step-by-step  integration  of  the 
turbulent  boundary  layer;  K  is  1  in  tlie 
predictor  stage  and  2  in  the  corrector 
stage 
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FORTRAN 

Variable 

Name 

Variable 

Definition 

L 

L 

Control  variable,  meaningful  only  when 
the  transition  criterion  of  Reference  2  is 
being  used  with  J  •  1;  L  =  1  gives  a  flat 
line  for  RTHT-RTHN  at  large  values  of  TP 
and  L  1  gives  a  sloping  line  for  it 

NCPTS 

Number  of  the  point  downstream  of  which 
the  point  spacing  begins  to  increase 

NCI 

NCPTS  +  1 

NC2 

NCPTS  +  2 

NC3 

NCPTS  +  3 

NC4 

NCPTS  +  4 

NMl 

NETS  -  1 

NM2 

NPTS  -  2 

NM3 

NETS  -  3 

NM4 

NPTS  -  4 

NM5 

NPTS  -  5 

NM6 

NPTS  -  6 

NETS 

N 

Total  number  of  points 

NPl 

NPTS  +  1 

0ML2 

9.l\} 

Momentum  area  divided  by  square  of  total 
length 

0ML2D 

OML2  far  downstream  of  the  body 

0ML2E 

OML2^®^ 

0ML2  at  tail 

0ML2T 

OML2^^^ 

OML2  at  the  transition  point;  a  function 
of  PGP  used  to  find  the  neutral-stability 

point 


101 


FORTRAN 

Variable 

Name 


Variable 


Def Inltion 


PGP 

.  ,  ,  .  RTH2RL  *  DUDX 

A  function,  equal  to  *  „  *  u  • 

used  to  find  the  neutral  stability  point 
and,  if  it  occurs,  laminar  separation; 

2 

PGP  is  also  equal  to - where  0,  V 

^  V  ds’ 

4  .•  4  4 

and  —  are  dimensional 
ds 

PGPN 

PGP^"^ 

PGP  at  the  neutral-stability  point 

PI 

7T 

Ratio  of  circumference  of  a  circle  to 
its  diameter 

PSI12 

Entrainment  area  divided  by  square  of 
length 

PSIL2E 

PSIL2^®^ 

PSIL2  at  tail 

RL 

Reynolds  number  based  on  total  length  L 
and  U^,  the  velocity  far  upstream 

RS 

R 

s 

Reynolds  number  based  on  distances  along 
the  body  and  U^,  the  velocity  far  upstream 

RTH 

"e 

Reynolds  number  based  on  momentum  thick¬ 
ness  and  local  velocity  outside  the 
boundary  layer;  RTH  is  equal  to  both 

0IJ 

— ,  where  9,  U,  and  V  are  dimensional, 

and  to  TH  *  U  *  RL 

RTHN 

RTH  at  neutral-stability  point 

RTHT 

RTH  at  transition  point 

RTH2RL 

RTH^/Rj^ 

SEGA 

sec  a  (1) 

Secant  of  ot,  evaluated  at  I;  a  is  the 
angle  between  the  tangent  to  the  body 
surface  and  the  x-axis 
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FORTRAN 

Variable 

Name 


Variable 


Definition 


SECl 

SEGA  (I  +  1) 

Secant  of  a,  evaluated  at  I  +  1;  a  is 
the  angle  between  the  tangent  to  the 
body  surface  and  the  x-axls 

SC 

0 

SG  =  /2/C  ,  where  C  is  the  local  shear- 
1  T 

stress  coefficient 

(e) 

ELSL2E 

SHE 

SH^ 

SH  at  tail,  SHE  -  q^l2E 

SKIP 

SKIP 

Variable  used  to  control  point  at  which 
reading  of  a  tape  begins 

TF 

TF 

A  function  of  TP,  used  in  the  transition 
criterion  of  Reference  2,  according  to 
which  transition  occurs  at  that  point 
downstream  of  neutral  stability  where 

RTH  -  RTHN  first  exceeds  TF;  there  are 
three  alternative  functional  relationships 
between  TF  and  TP,  depending  on  the  values 
of  J  and  L 

TFR 

TFR 

TF  +  RTHN 

TFT 

TF^t) 

TF  at  transition  point 

TF53 

Tf(53) 

2  *  ELD  *  (Y(I)  *  RTH  -  YN  *  RTHN); 
according  to  the  transition  criterion  of 
Reference  1,  transition  occurs  where 

TF53  first  exceeds  TP53 

TH 

9/L 

Momentum  thickness  divided  by  overall 
length 

THE 

TH^e) 

TH  at  tail 

THT 

TH^t) 

TH  at  transition  point 

TP 

TP 

A  function  of  body  geometry,  defined  as 

TP  =  L  ,  dY(I) 

Y(I)  *  ELD  dX 

TPT 

Tp(t) 

TP  at  transition  point 
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FORTRAN 

Variable 

Variable 

Def inition 

Name 

TP53 

Tp(53) 

Fmictioti  used  lo  predict  transition 
described  in  Kiderence  1  and  defined  as 

TP53  =  450  +  4()0e*^'''^ 

U 

Velocity  in  flow  field  immediately  out¬ 
side  bound.'irv  I.ivit,  divided  by  velocity 
far  upstn-.im;  dimensioned  to  be  a  func¬ 
tion  t)f  I 

UE 

,(o) 

U  at  tail 

UN 

,(n) 

U  at  nout ra 1 -s tab i 1 ity  point 

UT 

u(t) 

U  at  transition  point 

UU 

UU 

U  evaluated  at  center  of  integration  step 
by  linear  interpolation;  dimensioned  to 
be  a  function  of  1 

VL3 

V/L^ 

Body  volume  divided  by  cube  of  total  length 

X 

x/L 

Axial  coordinate  divided  by  total  length; 
dimensioned  to  he  a  function  of  I 

XN 

x<"> 

X  at  neutral-stability  point 

XT 

X  at  transition  point 

XTRIP 

j^(trip) 

If  transition  has  not  occurred  upstream 
of  XTRIP  and  if  ITRIP  is  set  equal  to  1, 
transition  is  forced  to  occur  at  XTRIP 

XX 

Axial  coordinate,  identical  to  X  except 
that  it  is  not  dimensioned  as  a  function 
of  I 

Y 

Body  radius  divided  by  total  length; 
dimensioned  tc'  be  a  function  of  I 

YN 

Y(n) 

Y  at  neutral-stability  point 

YT 

yCt) 

Y  at  transition  point 

104 


LISTING  AND  SAMPLE  OUTPUT  OF  DPOUT 
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